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ABSTRACT 


This  thesis  summarizes  an  investigation  on  the  effect  of  an  ice 
cover  on  open  channel  surges.  The  effect  of  ice  on  total  channel  friction 
is  analyzed  first  and  this  is  then  applied  to  open  channel  surges. 

The  propagation  of  surges  is  calculated  by  means  of  the  method 
of  characteristics.  With  the  use  of  numerical  techniques,  the  computer 
is  employed  to  solve  the  characteristic  equations.  The  behavior  of  surges 
is  tested  under  various  channel  conditions. 

A  radical  change  in  the  profile  and  propagation  velocity  of  a 
surge  is  observed  with  the  addition  of  an  ice  cover  to  a  channel. 
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GLOSSARY  OF  SYMBOLS 


cross-sectional  area  of  flow. 

top  width  of  flow. 

wave  propagation  velocity. 

positive  characteristic  curve. 

negative  characteristic  curve. 

depth  of  flow  under  open  water  conditions. 

depth  of  flow  under  ice  covered  conditions. 

friction  factor. 

water  surface  fall. 

acceleration  due  to  gravity. 

stage. 

sliding  coefficient  -  ratio  of  the  discharge  under  ice 
to  open  water  discharge  at  the  same  stage. 

relative  roughness  height. 

points  of  known  conditions. 

Manning's  roughness  factor  for  the  total  cross-section. 
Manning's  roughness  factor  for  ice  boundary. 

Manning's  roughness  factor  for  bed  boundary, 
point  of  unknown  conditions, 
wetted  perimeter. 

•  i 

d i scharge 

hydrau 1 i c  rad i us  . 

water  surface  slope. 

slope  of  energy  grade  line. 


bed  slope. 
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time. 


bouyant  thickness  of  ice. 
velocity  at  a  point, 
shear  velocity. 

average  cross-sectional  velocity, 
distance  along  channel. 

measurement  of  depth  at  point  velocities. 

distance  from  ice  boundary  to  point  of  maximum  velocity 

distance  from  bed  boundary  to  point  of  maximum  velocity 

channel  bottom  elevation. 

constants. 

unit  weight  of  water. 


bed  shear  stress. 


1. 


CHAPTER  I 

INTRODUCTION 

1 .  1  The  Problem 

On  northern  rivers,  dams  are  required  for  the  regulation  of  flow 
and  the  generation  of  hydro-electric  power.  To  meet  these  demands,  which 
are  greatest  during  the  winter,  the  release  of  large  quantities  of  water 
is  required.  The  hydro-electric  power  is  used  for  the  peak  daily  electrical 
loads  thus,  resulting  in  the  periodic  release  of  water.  This  rapid  change 
in  the  flow  rate  causes  a  surge  to  form  and  propagate  through  the  channel 
downstream  from  the  dam.  Since  the  advent  of  the  high-speed  digital  com¬ 
puter,  there  has  been  considerable  use  of  numerical  techniques  for  the 
solution  of  engineering  problems.  Recently  these  methods  have  been  used 
in  the  simulation  of  open  channel  wave  problems.  However,  they  have  only 
been  applied  to  short  channels  having  a  free  water  surface  condition.  In 
areas  of  cold  winters  the  rivers  and  streams  are  covered  with  ice  having 
appreciable  thickness  and  structural  strength.  These  numerical  techniques 
have  not  as  yet  been  applied  to  surges  under  an  ice  cover. 

The  present  investigation  was  carried  out  in  an  attempt  to  develop 
a  means  of  predicting  the  behaviour  of  surges  under  these  conditions. 
Presently,  in  Alberta,  the  releases  from  the  Brazeau  Dam  propagating  down 
the  Brazeau  and  North  Saskatchewan  Rivers  are  of  particular  interest.  More 
dams  are  planned  for  these  areas  resulting  in  a  greater  need  for  a  knowledge 
of  the  behavior  of  surges  under  ice  cover. 
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1 . 2  Purpose  and  Scope 

The  purpose  herein  was  to  simulate,  on  the  digital  computer,  surges 
propagating  down  an  ice  covered  prismatic  channel.  From  this  simulation 
the  effect  of  a  surface  ice  cover,  having  structural  strength,  on  the 
behavior  of  open  channel  surges  could  be  studied.  The  main  objectives 
of  the  study  were: 

1.  To  determine  the  effect  of  channel  roughness  on  open 
channel  surges. 

2.  To  determine  the  physical  behavior  of  open  channel  surges 
under  surface  ice  cover  conditions. 

In  order  to  simulate  an  ice  covered  channel  it  was  necessary  to 
determine  the  effect  of  the  ice  on  the  roughness  of  the  channel.  There¬ 
fore,  the  investigation  initially  dealt  with  ice  roughness  and  later 
surges  were  considered.  The  scope  of  the  investigation  was  limited 
by  the  lack  of  time  and  equipment  for  obtaining  roughness  measurements  as 
well  as  the  limit  on  computer  time  in  order  that  the  simulation  be  econom¬ 
ically  feasible  for  future  use. 

Data  on  the  channel  roughness  was  the  first  to  be  obtained.  Con¬ 
tinuous  stage  readings  were  recorded  throughout  the  ice  cover  period. 

These  were  related  to  discharge  readings  taken  at  a  section  of  the  river 
which  was  ice  free.  During  mid-winter  velocity  profiles  were  taken  to 
obtain  additional  data  on  the  effects  of  ice  cover  on  the  total  roughness 
of  the  channel . 

Using  numerical  techniques  and  the  digital  computer,  a  simulation 
of  surges  propagating  down  a  channel  was  made.  Runs  were  carried  out  under 
different  flow  rates  and  different  channel  roughnesses.  As  many  as 
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possible  of  these  runs  were  compared  to  actual  data  available  on  surges  on 
the  North  Saskatchewan  River. 


CHAPTER  I  I 


FRICTION  FACTOR  FOR  STEADY  FLOW 
UNDER  ICE  COVER 

2 .  1  Genera  I 

For  a  given  discharge,  an  ice-covered  river  maintains  a  higher 
stage  than  an  uncovered  river.  Part  of  this  increase  in  stage  is  due  to 
the  bouyant  thickness  "t."  of  the  ice  cover.  The  free  surface  of  the 
water  is  raised  approximately  0.9t.;  the  numerical  coefficient  depends  on 
the  ice  density.  An  additional  stage  increase  is  caused  by  the  increase 
in  frictional  resistance  due  to  the  roughness  of  the  underside  of  the  ice. 

The  effect  of  the  increase  in  friction  is  to  reduce  the  velocity  from  that 
of  the  uncovered  channel  wi th  an  equivalent  stage.  To  compensate  for  the 
decrease  in  velocity,  the  cross-sectional  area  must  increase,  thus  increas¬ 
ing  the  stage.  The  total  increase  in  stage  due  to  the  ice  cover,  from  the 
stage  for  the  same  open  water  discharge,  is  usually  referred  to  as  back¬ 
water.  This  backwater  changes  during  the  ice  cover  period  since  the  ice 
is  a  variable  feature  in  the  channel.  During  freeze-up,  there  is  a  large 
but  undefined  increase  in  friction.  This  friction  decreases  and  is  rather 
stable  during  midwinter.  Preceeding  break-up  there  is  again  a  poorly 
defined  rise  in  frictional  effects.  Some  methods  of  relating  the  stage 
to  discharge  have  been  proposed  but  all  give  limited  success. 

2 . 2  Literature  Review 

A  review  of  the  literature  indicates  that  very  little  work  has 
been  done  in  this  field  due  to  the  many  variables  involved  and  the  difficulty 
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in  obtaining  data  during  ice  cover.  The  majority  of  the  work  only  considers 
a  general  increase  in  stage  due  to  the  ice  cover,  with  little  consideration 
given  to  the  causes  of  this  increase. 

Tillinghast  (1905)  did  a  study  on  six  rivers  which  had  varying  dis¬ 
charges  from  100  to  6000  cfs .  From  his  studies,  he  found  that  it  was  more 
reliable  to  use  gauge  to  the  water  surface  rather  than  to  the  under  surface 
of  the  ice.  The  sliding  coefficient  "k",  which  is  the  ratio  of  measured 
discharge  to  open  water  discharge  at  the  same  stage,  was  found  to  vary 
between  0.30  and  0.75  throughout  the  season. 

In  studying  the  complex  manner  in  which  ice  forms,  Hoyt  (1913)  found 
that  each  stream  presented  different  conditions,  therefore  making  it  diffi¬ 
cult  to  formulate  one  method  for  calculating  discharge.  He  proposed  three 
methods : 

(i)  Water  surface  gauge  readings  may  be  applied  directly  to  the 
open  water  rating  curve  provided  the  stream  is  open  at  the 
control  section  and  no  backwater  exists  at  the  gauge. 

(ii)  Water  surface  gauge  readings  may  be  applied  directly  to  a 

special  rating  curve  based  on  winter  discharge  measurements, 

(iii)  Discharge  measurements  may  be  used  in  connection  with  gauge 

heights  and  data  showing  climatic  conditions  and  the  occurrence 
of  ice. 

In  using  a  backwater  correction,  Kolupaila  (1938)  saw  a  correlation 
between  the  correction  and  air  temperature.  He  concluded  that  "the  rela¬ 
tion  between  winter  and  summer  discharges  depends  on  three  factors:  (i) 
the  pecul ar i ties  of  the  roughness  of  the  ice  and  bed  (ii)  the  different 
slopes  of  the  water  level  and  (iii)  the  thickness  of  the  ice  in  relation 
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to  the  depth."  He  found  that  k  drops  from  1.0  to  0.2  or  0.3  during  freeze- 
up,  then  increases  slowly  all  winter  and  rapidly  rises  to  1.0  as  break-up 
occurs.  He  proposed  that  three  measurements  of  discharge  be  made,  the 
first  right  after  freeze-up  and  two  subsequent  ones  during  the  winter, 
thus  enabling  the  value  of  the  k  to  be  determined. 

Anchor  ice  is  ice  which  forms  on  underwater  rocks.  This  formation 
is  a  result  of  cooling  the  underwater  surfaces  below  the  freezing  point  by 
a  slight  supercooling  of  the  water.  The  coating  of  anchor  ice  may  be 
several  inches  thick  thus  causing  the  stage  to  rise.  An  increase  in  stage 
due  to  anchor  ice  will  not  occur  in  ice  covered  channels  since  the  ice 
cover  effectively  cuts  down  the  supercooling  of  the  water  (Chow  196A). 

Carey  ( 1967)  observed  ripplelike  and  dunelike  features  on  the  under¬ 
side  of  the  ice  cover.  From  his  observations  he  developed  two  methods  of 
calculating  discharge:  the  s tage-f a  1 1  - rel at i on  method  and  the  pipe-flow 
method.  The  first  method  basically  involves  two  graphical  relationships: 

(1)  a  relation  between  stage  and  discharge  for  some  fixed  condition  of  fall 
of  the  water  surface  and  (2)  a  relation  between  discharge  ratios  and  corres 
ponding  fall  ratios.  If  the  fall  is  not  constant  with  stage  an  additional 
graphical  relationship  is  required  between  stage  and  fall.  The  method 
depends  on  the  principle  that: 

=  function 

h 

where  Q  and  Q2  are  different  discharges  and  F]  and  F^  are  the  correspond¬ 
ing  falls  of  the  water  surface.  The  pipe-flow  method  makes  use  of  a  modi¬ 
fied  Darcy-Wei sbach  equation: 


8g 


f  mod 
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where  f  mod  is  a  modified  friction  factor,  S  is  the  water  surface  slope 
and  AR  is  the  section  factor.  A  modified  friction  factor  is  used  to 
compensate  for  the  fact  that  the  water  surface  slope  was  used  rather  than 
the  energy  gradient.  From  the  graphical  relationships  of  modified  friction 
factor  versus  water  surface  slope  and  section  factor  versus  gauge  height, 
the  discharge  is  obtained.  The  stage-fall  method  is  the  more  attractive 
of  the  two  because  it  is  more  accurate  and  it  requires  less  field  and 
of f i ce  work. 

Some  studies  have  been  done  on  the  calculation  of  a  composite 
roughness  for  an  ice  covered  channel,  as  part  of  the  increase  in  the 
stage  is  due  to  the  increase  in  roughness.  In  1964,  Nezki khovski y  recom¬ 
mended  three  formulae  for  calculating  a  composite  roughness  factor  for 
natural  channels: 

Weighted  average  n  =  n]+n2  (23) 


Pavol ovski y 1 s  n 


(2.4) 


Belokon's  -  Sabaneer's  n  =  0.63n^ 

The  weighted  average  method  is  correct  only  in  form  as  it  does  not  relate 
the  roughness  to  the  hydraulic  radius.  Pavlovskiy  assumed  that  the  total 
hydraulic  radius  was  equal  to  the  mean  depth  and  also  the  hydraulic  radii 
of  the  two  different  parts  of  roughness.  He  recognized  his  shortcomings 
and  wrote  that  his  formula  is  "useful  only  as  a  first  approximation". 

The  Belokon  -  Sabaneer  method  is  the  most  accurate.  It  is,  however,  de 
fective  when  the  ice  cover  is  ideally  smooth  and  therefore  it  must  only 
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be  applied  to  closed  channels  with  some  roughness.  With  the  above  formu¬ 
lae,  the  discrepancy  in  the  results  of  the  computation  of  n  with  observed 
values  of  n^  and  n^  was  less  than  S%. 

In  his  studies  of  two  power  channels  in  Sweden,  Larsen  (1969) 
developed  a  composite  roughness  formula: 

(y,/y,)5/3  I  +  i 

J_=  n  j  n2  -------------  -  (2.6) 

"  (l/2)2/3  (y,/y2+.)5/3 

where  and  y^  are  the  distances  from  the  bottom  of  the  ice  to  the  point 

of  maximum  velocity  and  from  the  bed  to  the  point  of  maximum  velocity, 

respectively.  Assuming  a  logarithmic  velocity  distribution,  he  determined 

a  roughness  coefficient  "k  11  for  each  part  of  the  boundary.  The  roughness 

*  1  /6 

coefficient  was  then  related  to  Manning's  n  by:  n  =  ck  .He  also 

® 

observed  bed  forms  on  the  underside  of  the  ice  cover. 

Recent  studies  have  been  carried  out  on  the  Peace  River  and  Slave 

River  in  northern  Alberta  by  the  Arctic  River  Work  Group  of  the  Federal 

Department  of  Energy,  Mines  and  Resources  in  Calgary.  In  1966,  Campbell 

realized  that  a  single  rating  curve  could  not  be  used  for  the  full  duration 

of  a  winter.  A  two  part  rating  curve  (named  the  W^  and  W^  method)  was 

developed,  with  the  periods  chosen  being,  "freeze-up  to  the  lowest  winter 

gauge  height"  and  "lowest  winter  gauge  height  to  break-up".  In  comparing 

six  methods,  Morton  (1967)  found  that  the  "Normal  Backwater"  method,  which 

uses  a  backwater  correction  to  the  open  channel  rating  curve,  gave  the 

best  accuracy  for  winter  streamflow.  Bennett  (1968)  modified  the  and 

W  method  to  give  the  discharge  equation: 

B 

2 

Q  =  Aq  +  Aj  (stage)  +A2(stage)  +  A^  (coded  date)  +  A^ (cumu 1  at i ve 
degree  -  days  since  freeze-up)  ------------  -(2.7) 
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He  suggested  that  three  discharge  measurements  be  made  during  the  winter 
to  determine  the  values  of  A:  (i)  immediately  after  freeze-up  (ii) 
approximately  the  time  of  minimum  gauge  height  (iii)  close  to  spring 
b  reak-up . 

Thus  far,  the  major  portion  of  the  ice  studies  have  been  done 
during  the  middle  portion  of  the  winter.  There  has  been  some  considera- 
tion  of  the  break-up  period  but  the  freeze-up  period  has  been  omitted 
since  the  backwater  is  generally  very  erratic  during  this  time.  Climatic 
conditions,  large  increases  in  friction  and  ice  jamming  are  all  causes 
of  the  erratic  backwater.  It  is  unfortunate  that  this  period  has  been 
omitted  since  its  effects  contribute  to  the  condition  of  flow  during  the 
remainder  of  the  ice  cover  period. 

2 . 3  Theoretical  Aspects  of  Determining  Roughness  from  Measured 

Velocity  Distribution. 

The  forming  of  an  ice  cover  on  a  wide  open  channel  causes  a  radical 
change  in  its  hydraulic  conditions.  As  a  result  of  the  ice  cover,  the 
open  channel  becomes  a  closed  conduit.  Its  cross-sectional  area  adjusts 
to  the  need  since  the  ice  is  a  floating  mass.  This  surface  cover  produces 
a  large  increase  in  the  wetted  perimeter.  Generally  the  roughnesses  of 
the  two  boundaries  are  not  the  same.  Therefore,  it  becomes  necessary  to 
compute  a  composite  roughness  to  determine  the  head  losses. 

Assuming  the  Prandtl -Karman  equation  for  fully  developed  turbulent 
flow  to  be  valid,  the  roughness  of  the  boundary  can  be  determined  from 
the  velocity  distribution.  The  velocity  distribution  near  the  boundary 
is  governed  by  the  roughness  and  distance  from  the  boundary  according  to: 
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U: 


a  lo9  y/k  +  3 


(2.8) 


From  the  experimental  work  of  Nikuradse  the  values  of  a  and  3  are  5.75 
and  8.5  respectively.  Plotting  the  measured  velocities  on  semi - 1 ogar i th- 
mi c  paper  should  produce  a  straight  line.  The  slope  of  the  line  is  equal 
to  5.75*  U* .  Choosing  U  and  a  corresponding  y,  the  equivalent  roughness, 
k^,  can  be  determined  by  substitution  into  equation  (2.8). 

For  rough  turbulent  flow  in  a  wide  channel,  Keulegan  (1938)  showed 

that: 

IT  =  6.0  +  5.75  log  __y  - - - - - (2.9) 

U*  k 

s 

Hanning's  formula  for  a  wide  channel  can  be  written  as: 

U  =  1.486  y2/3  So,/2  .  (2.10) 

n 

with  U*  equal  to  /  g7So  equation  (2.10)  can  be  rewritten  as: 

U_  =  1,486  1  y,/6  . . (2.11) 

U*  n  77“ 

Equating  equations  (2.9)  and  (2.11)  results  in: 

n  =  _ 1.486  y1/6  - . (2.12) 

/g  (6.0  +  5*75  log  y /,  ) 

Ks 

From  the  calculated  values  of  y  and  ks>  a  value  for  Manning's  n 


can  be  calculated. 


*■  '^j 

- 


CHAPTER  I  I  I 


ANALYSIS  OF  ROUGHNESS  UNDER  ICE  COVER 

3  •  1  General 

A  preliminary  study  of  the  ice  roughness  on  the  North  Saskatchewan 
River  through  Edmonton  was  carried  out  during  the  winter  of  1969  -  70. 

The  data  collected  consisted  of  continuous  stage  recordings  taken  through¬ 
out  the  ice  cover  period  and  velocity  profiles  taken  during  mid-winter. 
From  this  information  the  variation  in  roughness  during  the  winter  can  be 
readily  detected.  A  roughness  coefficient  was  calculated  for  the  mid¬ 
winter  period.  The  data  collected  was  not  sufficient  to  give  a  complete 
analysis  of  the  additional  roughness  due  to  the  ice  cover,  but  it  gives 
a  fair  estimate  of  the  roughness  coefficient. 

3 •  2  Col  1 ect i on  of  Data 

At  three  sections  in  the  ice  covered  portion  of  the  river,  stage 
readings  were  taken  by  continuous  stage  recorders.  Two  of  the  recorders 
were  located  in  the  wells  of  pumping  stations  while  the  other  was  located 
in  a  bridge  pier.  Figure  3-1  gives  a  plan  of  the  North  Saskatchewan  River 
through  Edmonton  with  the  locations  of  the  recorders  labelled.  The  gauge 
located  in  the  pier  of  the  Quesnel  Bridge  operated  only  a  short  time  after 
freeze-up  and  then  it  became  frozen.  Continuous  pumping  at  varying  rates 
existed  at  the  California  Standard  Pumphouse.  The  unsteady  pumping  caused 
the  well  to  have  a  varying  amount  of  drawdown;  therefore,  inducing  error 
into  the  stage  recordings.  The  gauge  located  at  Mayfair  Park  gave  the 
most  accurate  readings  of  the  three  gauges.  The  occasional  pumping  which 
occurred  at  this  gauge  could  be  detected,  thus  the  stage  readings  with 
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FIGURE  3-1  MAP  OF  THE  NORTH  SASKATCHEWAN  RIVER  THROUGH  EDMONTON 
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drawdown  were  discarded. 

The  mean  daily  discharges  were  obtained  from  the  Water  Survey  of 
Canada  which  operates  a  continuous  stage  recorder  at  the  Low  Level  Bridge 
(See  Figure  3“1  for  the  location).  This  gauge  was  situated  in  a  section 
of  the  river  which  was  free  from  ice  cover.  Warm  water,  which  had  been 
used  as  cooling  water  for  the  power  plant,  discharged  into  the  river 
keeping  it  free  of  ice.  With  no  ice  cover,  the  discharges  were  simply 
obtained  from  the  open  water  rating  curve  (see  Appendix  A). 

The  velocity  distribution  in  a  cross-section  was  measured  at  the 
Mayfair  Park  gauge  during  mid-winter.  For  the  insertion  of  the  current 
meter,  six  inch  diameter  holes  were  drilled  through  the  ice.  Velocity 
distributions  were  taken  every  50  feet  with  some  additional  measurements 
taken  near  each  bank  (see  Appendix  A).  A  small  Ott  current  meter,  which 
included  a  counter  and  timer,  was  used  for  taking  the  velocity  measurements. 
At  each  hole,  velocity  readings  were  taken  at  every  0.25  feet  for  the 
entire  depth.  In  addition  to  the  velocity  measurements  at  each  hole;  the 
river  depth,  ice  thickness  and  the  water  level  with  respect  to  the  ice 
surface  were  noted.  The  above  measurements,  with  the  exception  of  velocity 
readings,  were  also  made  at  the  California  Standard  Pumphouse. 

3-3  Treatment  of  Stage-Discharge  Data 

The  ice  cover  period  on  the  North  Saskatchewan  River  lasted 
approximately  four  months;  from  November  21,  1969  to  March  21,  1970. 

In  Figure  3-2,  a  stage-discharge  relation  was  plotted  for  the  Mayfair 
Park  gauge.  Rating  curves  for  three  different  ice  cover  periods  were 
obtained:  (i)  freeze-up  (ii)  break-up  (iii)  mid-winter.  The  open 

water  curve  is  also  shown  on  this  plot. 
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FIGURE  3-2  STAGE- DISCHARGE  CURVES 
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Of  the  three  periods,  the  freeze-up  period  requires  the  highest 
stage  for  a  given  discharge.  The  freeze-up  period  was  preceeded  by  an 
ice  jam  just  below  the  Low  Level  Bridge.  The  ice  jam  induced  additional 
errors  to  the  discharge  measurements  by  causing  increases  in  the  stage 
due  to  ice  debris  floating  in  the  river,  a  backwater  from  the  ice  jam 
and  a  decrease  in  the  width  of  the  river.  As  time  progressed  during 
the  freeze-up  period,  the  stage  for  a  given  discharge  rapidly  decreases. 
This  was  due  to  the  "washing  out"  of  the  ice  jam  and  the  smoothening  of 
the  underside  of  the  ice.  The  ice  thickness  increases  rapidly  during 
this  period  but  it  can  be  seen  from  the  plot  that  the  decrease  in  stage 
due  to  the  smoothening  of  the  ice  more  than  offsets  the  increase  in 
stage  due  to  the  increased  ice  thickness. 

The  mid-winter  stage-discharge  relationship  is  stable.  During 
this  period  the  ice  grows  very  little  in  thickness  and  the  decrease  in 
roughness  is  small.  This  period  continues  for  the  major  portion  of  the 
ice  cover  period  as  it  lasts  for  all  but  the  approximate  three  week 
period  of  freeze-up  and  one  week  of  break-up. 

The  break-up  period  sees  a  rise  in  the  stage  for  a  given  discharge 
During  this  period  the  thinnest  ice  sections  of  the  river  and  the  small 
tributary  streams  break  up.  This  broken  ice  then  flows  under  the  main 
sheet  of  ice,  thus  increasing  the  effective  roughness.  The  additional 
frictional  losses  result  in  the  increased  stage. 

The  stage-discharge  curve  for  an  ice  covered  channel  is  a  variable 
feature.  The  changes  in  roughness  and  ice  thickness  cause  this  variation 
Ice  jams  do  not  occur  as  regular  as  the  above  causes  but  thei r  effect 
can  be  very  large  due  to  changes  in  channel  width  and  backwater  effects. 
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The  length,  of  the  three  periods  previously  discussed,  varys  from  year 
to  year  depending  on  the  climatic  conditions.  A  measurement  must  be 
made  during  each  period  to  determine  the  stage-discharge  relationship 
for  ice  cover. 

3 •  4  Treatment  of  Velocity  Data. 

The  velocity  distribution  was  measured  at  the  Mayfair  Park  gauge 
on  February  7,  1970.  The  ice  was  approximately  2  1/2  months  old  with 
an  average  thickness  of  1  1/2  feet.  The  measured  velocity  distribution 
is  shown  in  Figure  3“3. 

The  velocity  distribution  measurements  are  treated  according  to 
the  theoretical  approach  using  the  Prandtl  -  Karman  law.  Measurements 
were  made  on  12  verticals  and  a  regression  analysis  was  utilized  to 
obtain  the  best  fit  straight  line  through  the  points.  The  equivalent 
roughness  height,  kg,  was  then  found  and  related  to  Manning's  n  according 
to  equation  (2.12).  Once  the  Manning  coefficient  was  obtained  for  both 
the  bed  and  the  ice  boundaries,  a  composite  value  was  calculated  by  the 
methods  mentioned  in  the  previous  chapter.  The  calculated  data  are  tab¬ 
ulated  and  explained  in  Appendix  B.  The  values  obtained  for  the  composite 


n  are: 

Larsen - 0.0219 

Weighted  Average  -----------------  0.0225 

Belokan  and  Sabaneer  ---------------  0.0231 

Pavlovskiy  --------------------  0.0238 


The  ice  cover  was  found  to  be  smooth  with  a  Manning  coefficient 
of  0.0150.  From  observations  made,  the  ice  taken  from  the  channel 
appeared  to  be  completely  smooth.  This  is  contrary  to  what  was  observed 
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FIGURE  3-3  CROSS-SECTION  OF  THE  NORTH  SASKATCHEWAN  RIVER 
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by  both  Carey  (1967)  and  Larsen  (1969).  The  ripplelike  configurations 
observed  by  Carey  were  very  small  and  they  did  not  become  noticeable 
until  late  in  the  winter.  The  author  believes  that  these  ripplelike 
and  dunelike  features  may  be  caused  by  a  phenomenon  similar  to  that 
of  sand  dunes.  The  observations  made  on  the  North  Saskatchewan  were  in 
the  plain  bed  region  while  those  made  by  Larsen  were  in  the  dune  phase 
due  to  greater  depths. 

It  should  be  noted  that  the  values  of  equivalent  roughness  of 
the  bed  showed  rather  large  fluctuations.  This  would  be  expected  since 
the  bed  was  mostly  small  gravel  with  some  rather  large  rocks  spread  on 
top.  The  roughness  of  the  ice  was  considerably  more  regular. 

The  flow  depth  in  an  ice  covered  channel  is  greatly  increased 
over  that  of  an  open  channel  due  to  the  increase  in  the  wetted  perimeter 
The  percentage  increase  will  depend  on  the  roughness  of  the  ice  and  the 
original  roughness  of  the  channel. 


CHAPTER  IV 


SURGES 

k .  1  Genera  1 

This  chapter  briefly  discusses  the  previous  work  done  on  surges 
as  well  as  the  theoretical  aspects  of  the  problem.  Open  channel  surges 
are  considered  first  and  then  modifications  are  made  to  include  the  ice 
cover  condition.  The  theory  is  so  developed  that  it  can  be  easily 
adapted  to  the  computer. 

A . 2  Literature  Review 

In  reviewing  the  literature  there  has  recently  been  a  considerable 
increase  of  interest  in  surges.  This  is  a  result  of  the  advent  of  the 
digital  computer  which  has  enabled  these  problems  to  be  solved  more 
rapidly.  Most  of  the  studies  have  been  concerned  with  changes  in  the 
flow  rate  and  tidal  effects  on  open  channels.  The  effect  on  surges 
caused  by  an  ice  cover  on  open  channel  has  only  been  briefly  mentioned. 

Grushevskiy,  ( 1 965 )  discusses  the  characteristics  and  peculiarities 
of  the  various  types  of  unsteady  movement  of  water  in  unerodable  natural 
channels.  He  describes  the  propagation  of  surges  under  varying  channel 
conditions.  In  mentioning  ice  cover  he  states  that  the  celerity  and 
height  of  surges  are  decreased  under  these  conditions.  There  are  no 
specific  conclusions  reached  in  his  paper  as  it  is  a  general  description 
of  unsteady  flow. 

Cooper  (1966)  performed  some  laboratory  experiments  on  the  effect 
of  a  surface  cover  on  open  channel  surges.  Using  surface  covers  of  vary- 
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ing  rigidity,  roughness  and  weight,  he  observed  surges  propagating  through 
initially  still  water.  From  his  experiments  he  found  that  the  addition  of 
a  surface  cover  reduces  the  velocity  and  height  of  surges.  The  decreases 
he  observed  in  laboratory  models  were  not  as  great  as  those  in  the  field. 

To  compute  surge  propagation  against  an  adverse  current,  Collins 
and  Fersht  (1968)  developed  a  numerical  scheme  for  the  integration  of  the 
basic  non-linear  equation  with  friction.  They  used  a  mixed  technique 
involving  finite  difference  forms  combined  with  a  fourth-order  Runge- 
Kutta  method.  They  checked  their  procedure  against  observations  made  dur¬ 
ing  a  hurricane.  The  computed  results  give  a  reasonably  good  representation 
of  the  actual  observations. 

In  1968,  Baltzer  and  Chintu  Lai  simulated  unsteady  flows  in  rivers 
using  three  different  methods.  They  used  the  method  of  characteristics, 
the  power  series  method  and  the  implicit  finite  difference  method,  all 
of  which  are  suited  for  computer  execution.  Comparisons  of  simulated 
flows  obtained  by  using  each  of  these  methods  with  numerous  field  measure¬ 
ments  indicate  generally  good  agreement. 

In  demonstrating  the  use  of  the  digital  computer  for  the  simulation 
of  prismatic  open-channel  wave  problems,  Martin  and  De  Fazio  ( 1 969)  used 
the  finite-difference  technique.  Varying  boundary  conditions  were  imposed 
on  the  channels  to  simulate  different  engineering  situations.  The  com¬ 
parison  of  the  computer  solutions  and  experimental  values  indicates  close 
agreement  especially  on  the  maximum  values  of  water-surface  elevation. 

A  tidal  study  was  carried  out  on  the  St.  Lawrence  River  in  1970 
by  Kamphuis.  He  also  applied  the  finite-difference  technique  to  solve 
the  equations  of  continuity  and  motion.  He  suggests  that  where  the  cross 
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sections  vary  considerably  with  distance  the  accuracy  of  the  solution  is 
much  improved  if  shorter  section  lengths  are  used.  The  computer  solutions 
were  similar  to  the  field  observations. 

The  previous  work  on  surges  has  mainly  considered  open  channel  surges. 
There  has  been  very  little  work  done  on  surges  under  an  ice  cover  and  it 
has  produced  only  limited  success.  The  computer  simulations  of  open  channel 
surges  have  provided  rather  good  results.  Although  several  conditions 
have  been  imposed  on  the  channels,  there  has  been  no  previous  attempt  to 
simulate  surges  under  an  ice  covered  channel. 

k . 3  Theoretical  Aspects  of  Surges 

Changing  the  rate  of  flow  causes  a  surge  to  form  and  propagate 
downstream  in  a  channel.  The  properties  of  the  surge  are  dependent  upon 
the  properties  of  the  flow  and  the  channel.  The  one-dimensional  equations 
of  continuity  and  linear  momentum  provide  the  necessary  tools  for  analyzing 
these  surges.  The  two  equations  allow  the  flow  characteristics  of  the 
channel  to  be  completely  expressed  in  terms  of  two  variables,  e.g.,  depth 
and  veloci ty . 

From  the  definition  sketch  in  Figure  4.1,  H  is  the  stage,  y  is  the 
depth  of  flow,  z  is  the  channel  bottom  elevation,  V  is  the  average  cross- 
sectional  velocity,  A  is  the  cross-sectional  area  of  flow,  B  is  the  top 
width  of  flow,  and  P  is  the  wetted  perimeter  of  the  open  channel.  If  x 
is  the  horizontal  co-ordinate  and  t,  time,  the  continuity  equation  may  be 
wr i tten  as : 

_1  CVA)  dx  =  'Bdx  lii  - ('|J) 

3x  3 1 


Simpl i fy i ng  gives : 
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FIGURE  4-1  DEFINITION  SKETCH 
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B  3H_  +  V  3A  +  A  3V_  =  0  - (4.2) 

3t  3x  3x 

A  3V_  +  3_H  +  V  3y  =  0  - (4.3) 

B  3x  3t  3x 


The  momentum  equation  is  obtained  by  equating  the  rate  of  increase 
of  momentum  in  the  control  volume  plus  the  rate  of  efflux  of  momentum  from 
the  control  volume  to  the  applied  force.  The  applied  force  is  a  sum  of 
the  pressure  force  and  the  bed  shear  force,  where  the  pressure  force  is 
-  yA  (3H/3x)  dx  and  the  shear  force  is  -  x  Pdx  =  -  yA  S^.  dx.  The  momen¬ 
tum  equation  is  thus  written  as, 


_3 
3 1 


1 

A 


y.  dxAV 
9 


_3_ 

3x 


y_  dxAV* V 


=  -yA  _3H_  dx 
3x 


3  (AV) 
3 1 


3  (AV-V) 
3x 


-g  3H_  -  g  S, 
3x 


yA  Sj-dx  -  -  (4.4) 
-  -  (4.5) 


The  left  hand  side  becomes 


A 


A  3 V 

3t 


+ 


V  3A  + 
3 1 


2AV  3V_  +  V2  3A 
3x  3x 


(h.6) 


Using  equation  (4.2)  and  multiplying  the  left  hand  side  by  A,  gives 


A  9V  +  V3A  +  2AV  3 V  +  V 
aT  3 1  3x 


-A  3V 
3x 


B  3H_ 
3 1 


- -  -(4.7) 


Subs  t i tut i ng: 

9A  =  3A  9H  =  B  3H_  - - - - - ^-8) 

3 1  3H  3 1  3t 

Expression  (4.7)  becomes 

A  3V  +  VB  3H  +  2AV  3V_  “  AV  3 V_  -  VB3H_  -  (4.9) 

ft  3T  3x  3x  3 1 
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which  reduces  to: 


A  9V  +  AV  3V 
3t  3x 


(4.10) 


Now  divide  express  ion (4. 1  0)  by  A  to  revert  back  to  the  left  hand  side  of 
equation  (4.5), 


aV_  +  V  8V_ 
8 1  8x 


(4.11) 


Thus  the  momentum  equation  reduces  to, 


_3V  +  V9V  +  g  3H  +  g  S  =  0 

at  3x  ax 


(4.12) 


Since  the  right  hand  sides  of  the  continuity  and  momentum  equations 
are  equal  to  zero  they  can  be  equated  by  using  X,  a  multiplying  constant: 


3V_+  V  9V  +  ga_H_+gS  +  X 
3 1  3  x  3  x 


A  av  + 

B  3x 


3H_ 

at 


v  ay 
ax 


=  0  -  (4.13) 


In  order  to  eliminate  H,  it  should  be  noted  that: 


H  =  z  +  y 


(4.14) 


3H_  =  3z_  +  3y_ 

ax  ax  ax 


(4.15) 


aH  =  ay 
at  at 


(4.16) 


Substituting  for  3H/ax  and  9H/3t  in  equation  (4.7)  gives 


av_  +  V  av_  -  gSo  +  g  3y  +  gSf  +  X 
at  ax  3x 


a  av_  + 
b  ax 


ay 

at 


V  ay 

ax 


-  -  (4 


where:  So  =  -az/3x 


(4.18) 
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Rearranging  equation  (4.17)  gives: 
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J<~  g(Sf-So)  -  -  -  (4.19) 


Now  find  X  to  make 

V  +  X  A  =  V  +  g/X  - - - - - (4.20) 

B 

which  is  X  =  +_  gB  -------------  -  (4.21) 

/  A 


The  terms  containing  9V/9t,  9V/9x,  9y/9t,  9y/9x  become  total  derivatives 


dx  =  V  +  X  A  =  V  +  g/X  - -  (4.22) 

dt  B 


dx  =  V  +  /  gA  '  =  V—  Co  - - -  - - (4.23) 

dt  /  B 


where  —  Co  =  +_  /  gA"  (4.24) 

/  B 


Therefore,  along  a  positive  characteristic  where: 

X  =  +  rji  and  dx_  =  V  +  Co  -------  -  (4.25) 

/A  dt 


equation  (4.19)  becomes: 

Co+  dV  +  dy  +  Co+  (Sf  -  So)  =  0  --------  (4.26) 

g  dz  dt 

Similarly  along  a  negative  characteristic  where: 

X  =  /  and  dx  =  V  -  Co  ----------  (4.27) 

/  A  dt 


equation  (4.19)  becomes: 

Co"  dV  -  dy  +  Co  (Sf  -  So)  =  0  ----------  (4.28) 

g  dt  dt 
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Figure  4-2  shows  the  characteristic  grid  on  the  x-t  plane  used  in 
solving  the  differential  equations.  With  subcritical  flow,  a  disturbance 
is  propagated  both  upstream  and  downstream  from  the  point  of  occurence, 
whereas  in  supercritical  flow  the  disturbances  only  propagate  downstream. 

This  can  be  seen  graphically  in  Figures  4~3  and  4-4  where  a  change  in  con¬ 
ditions  at  point  L  effects  all  the  shaded  area. 

4 . 4  Variations  For  Ice  Conditions 

The  hydraulic  characteristics  of  an  open  channel  are  radically  changed 
with  the  imposition  of  an  ice  cover.  The  flow  in  an  ice  covered  river  is 
in  some  ways  similar  to  the  flow  in  a  conduit.  In  a  wide  river  there  is 
one  important  difference  and  that  is,  that  the  cross-sectional  area  will  ad¬ 
just  according  to  the  need  because  the  ice  is  floating  on  the  water. 

The  ice  cover  produces  a  substantial  increase  in  the  wetted  perimeter, 
see  Figure  4-5.  Assuming  a  thin  cover  on  a  wide  channel,  the  hydraulic 
radius  is  reduced  to  near  half  the  value  of  the  open  channel  radius.  Also 
making  the  assumption  that  the  flow  is  uniform  and  keeping  the  flow  rate 
constant,  the  variation  in  depth  is  computed  from  Manning's  equation  as 
indicated  below.  For  an  open  channel 


(4.29) 


For  an  ice  covered  channel 


(,t.30). 


Equating  the  two  equations  gives: 
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FIGURE  4-2  GRID  OF  CHARACTERISTICS  IN  X-T  PLANE 
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FIGURE  4-3  CHARACTERISTICS  IN  SUBCRITICAL  FLOW 


FIGURE  4~4  CHARACTERISTICS  IN  SUPERCRITICAL  FLOW 
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HYDRAULIC  RADIUS 

OPEN  CHANNEL  R=A/P 
ICE  COVER  R  =  A/P*B 


FIGURE  4-5  HYDRAULIC  RADIUS  FOR  OPEN  AND  ICE 

COVERED  CHANNELS 
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If  the  ice  roughness  was  equal  to  the  bottom  roughness,  the  increase 
in  flow  depth  with  an  ice  cover  would  be  32%.  The  ice  is  generally  smoother 
than  the  bed,  therefore  a  more  practical  value  would  be  to  assume  the  ice 
roughness  coefficient  n  equal  to  half  the  bed  roughness  coefficient;  thus 
the  increase  in  depth  would  be  12%.  For  more  exact  results,  a  roughness 
based  on  both  boundaries  are  outlined  in  Chapter  II  should  be  applied  to 
determine  the  changes  in  depth  of  flow. 

The  above  modifications  were  applied  when  solving  the  characteristic 
equations  for  ice  cover  conditions.  Any  other  effects  due  to  the  ice  cover 
were  neglected. 


A . 5  Computer  Solution  to  Differential  Equations 

For  the  solution  of  the  differential  equations,  the  initial  conditions 
must  be  known  as  well  as  the  boundary  conditions  at  either  end  of  the  channel 
From  Figure  A-6,  the  conditions  of  time,  velocity,  depth,  and  distance  are 
known  at  points  L  and  M.  It  is  required  to  obtain  the  conditions  mentioned 
above  for  point  N.  Thus,  along  a  positive  characteristic,  equations  (A. 25) 
and  (A. 26)  become: 


tn  -  tl 

XN  '  XL 
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Co 


V  +  Co 

(\  -  vl\  .  -  Y' 


T  -  T 
N 


LJ 


N 


tn  -  tl 


(<1.33) 


+  Co+  (S  -  So)  =  0  - (4. 3*1) 


And  similarly  along  the  negative  characteristic,  equations  (4.27)  and  (4.28) 


become : 


K 
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FIGURE  4-6  DEFINITION  SKETCH 
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(4.36) 


The  values  used  for  Co  —  and  V  are  an  average  of  the  values  of  Co  — 

and  V  at  points  L  and  M.  The  four  unknowns  (Y..,  V..,  T,  Xkl)  at  point  N 

N  N  N  N 

can  be  obtained  from  the  four  equations.  This  is  a  quasi-linear  problem 
since  the  coefficients  Co  — ,  Co  ,  V  +  Co  and  V- Co  depend  on  the  values 
of  and  V  .  The  solution  to  these  equations  can  be  obtained  by  implement¬ 
ing  the  predictor  -  corrector  process. 


4. 6  Boundary  Conditions 


The  initial  condition  of  the  channel  is  imposed  on  the  x-axis  of 
the  x-t  diagram.  Generally  the  initial  condition  imposed  is  that  of  a 

steady-flow  condition,  although  any  condition  may  be  imposed  if  the 

velocities  and  depths  are  known  along  the  entire  x-axis. 

The  t-axis  at  one  end  of  the  channel  and  the  line  parallel  to 

the  t-axis  at  the  other  end  of  the  channel  are  the  lines  on  which  the 

boundary  conditions  are  imposed.  One  of  the  two  dependent  variables  must 
be  specified  as  a  boundary  condition.  At  the  boundaries  only  one  charact¬ 
eristic  leads  to  the  unknown  as  shown  in  Figures  4-7  and  4-8.  This 
results  in  two  differential  equations  available  for  the  solution  of  the 
condition  at  the  boundaries.  Since  either  depth  or  velocity  is  an  imposed 
condition  at  the  boundaries  and  the  distance  is  known,  the  two  available 
equations  provide  the  means  for  obtaining  the  two  unknowns. 
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FIGURE  4-8  RIGHT  HAND  BOUNDARY 
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Several  different  boundary  conditions  may  be  imposed  on  a  channel 
to  simulate  various  engineering  situations.  Some  of  the  more  typical 
boundary  conditions  are  hydrographs,  reservoirs,  dead  ends,  free-fall  and 
junctions.  With  the  initial  condition  and  boundary  conditions  imposed  on 
every  channel  that  comprises  a  total  system,  the  propagation  of  surges 
through  the  channels  can  be  simulated. 
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CHAPTER  V 

THE  SIMULATION  OF  SURGES 

5 . 1  Genera  1 

A  computer  program,  based  on  the  theory  expressed  in  the  previous 
chapter,  was  written  for  surge  propagation  through  a  channel.  This  pro¬ 
gram  was  tested  against  observations  made  on  the  Brazeau  and  North  Saskat¬ 
chewan  River  in  Alberta.  The  periodic  releases  of  large  quantities  of 
water  from  the  Brazeau  Reservoir  caused  surges  to  form  and  propagate  down 
the  channels.  Various  recording  stations  along  these  channels  provided 
information  on  the  surge  heights  and  speeds. 

Surges  of  various  sizes  and  frequency  were  imposed  at  the  upstream 
end  of  the  simulated  channel  and  were  allowed  to  propagate  through  the 
channel  under  various  conditions.  Computer  runs  were  made  on  both  open 
channel  and  ice  cover  conditions  with  the  solutions  being  provided  in  the 
form  of  plots.  The  velocities  and  depths  in  the  channel  were  plotted 
against  distance  at  various  times  and  against  time  at  various  stations 
a  1 ong  the  channel . 

5 . 2  Field  Observations 

In  order  to  test  the  accuracy  of  the  computer  solutions,  it  was 
necessary  to  obtain  some  field  data  on  surge  propagation.  The  required 
data  was  obtained  from  the  North  Saskatchewan  River  where  the  stage- 
recorders  on  the  river  record  the  propagation  of  the  surges,  which  are 
formed  from  the  controlled  releases  of  water  3t  the  Brazeau  Reservoir. 
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Figure  5  1  gives  a  plan  of  the  section  of  the  North  Saskatchewan 
and  Brazeau  Rivers  used  for  observations.  There  are  stage  recorders 
located  at  Rocky  Mountain  House,  Lodgepole,  Drayton  Valley  and  Edmonton 
on  the  North  Saskatchewan  River.  The  gauge  at  Drayton  Valley,  which  is 
operated  by  Calgary  Power,  gave  very  poor  readings  as  the  river  cross- 
section  was  not  stable;  this  resulted  in  the  discarding  of  those  readings. 
The  gauges  at  Edmonton  and  Rocky  Mountain  House,  which  are  operated  by  the 
Water  Survey  of  Canada,  record  the  stage  continuously  through  the  entire 
year.  The  river  at  the  Low  Level  Bridge  gauge  at  Edmonton  is  ice  free  for 
the  entire  year  whereas  an  ice  cover  forms  during  the  winter  on  the  river 
at  Rocky  Mountain  House.  The  cross-section,  recorder,  charts  and  rating 
curves  were  obtained  from  the  V/a ter  Survey  of  Canada  for  the  periods, 
September,  1967  to  March,  1968  and  May,  1969  to  May,  1970.  The  gauge  at 
Lodgepole,  which  only  records  during  the  summer  months,  was  a  recent  in¬ 
stallation;  therefore,  the  information  as  mentioned  above  was  available 
only  for  the  period  May,  1969  to  October,  1 969  - 

On  the  Brazeau  River  there  is  a  stage  recorder  located  just  below 
the  Big  Bend  Plant  on  a  section  of  the  river  which  is  ice  free  all  year. 
Since  recent  readings  were  not  accurate,  due  to  the  shifting  of  the  cross- 
section,  Calgary  Power  supplied  recorder  charts  and  rating  curves  for  the 
period  September,  1967  to  March, 1968.  The  data  obtained  from  this  gauge 
provided  very  good  information  on  the  variation  in  discharges  due  to  the 
releases  of  large  quantities  of  water  through  the  turbines  at  the  Big 
Bend  Plant. 

The  cross-sections  of  the  North  Saskatchewan  River  at  the  recording 
stations  were  found  to  increase  slightly  in  width  as  one  progressed  down- 
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stream,  but  they  could  be  approximated  by  trapezoidal  sections.  Similarly, 
from  the  longitudinal  profile  in  Figure  5"2,  it  can  be  seen  that  the  bed 
slope  slowly  decreases  as  one  progresses  downstream. 

5 . 3  The  Computer  Program 

Since  the  solution  of  surge  problems  are  readily  obtained  by  numeri¬ 
cal  techniques,  the  high-speed  digital  computer  is  very  well  suited  for  the 
task.  A  computer  program  was  written  in  Fortran  IV  language  using  the 
equations  developed  in  the  previous  chapter  with  a  separate  program  written 
for  plotting  the  results.  The  computer  analysis  was  carried  out  at  the 
University  of  Alberta  on  an  IBM  System/360,  Model  67,  computer  with  the 
results  being  plotted  by  a  model  663  Cal  Comp  Plotter. 

The  program  was  written  as  a  number  of  subroutines  so  as  to  increase 
its  flexibility  and  to  allow  for  the  ease  in  making  any  changes  which 
might  be  necessary.  The  tabulation  in  Table  5.1  gives  the  name  and  pur¬ 
pose  of  the  subroutines  with  a  complete  listing  of  program  available  in 
Appendix  C.  The  subroutines,  as  listed  in  the  appendix,  are  set  up  to 
handle  surges  in  an  ice  covered  channel.  Modifications  in  the  subroutines 
REACH,  BDY1  and  BDY2  were  required  in  changing  from  open  water  to  ice 
covered  channels.  The  changes  made  were  due  to  the  different  hydraulic 
radii  for  the  two  conditions.  This  was  previously  discussed  in  Chapter 
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FIGURE  5-2  LONGITUDINAL  PROFILE  OF  THE  NORTH  SASKATCHEWAN  RIVER 


TABLE  5-1  ~  NAME  AND  PURPOSE  OF  SUBROUTINES 


NAME 

PURPOSE 

WAV  21 

Control  Program 

SCRIB 

Writes  out  T,  X,  W,  Y  from  VECTR 

SORT 

Find  maximum  depth  and  maximum  and  minimum  velocities  for 
plot  scales. 

REACH 

Calculates  interior  points  on  X-T  diagram 

BDY  1 

Left  Boundary  Condition  (Discharge  versus  time) 

BDY  2 

Right  Boundary  Condition  (Stage-Discharge  Curve) 

TABL 

Interpolating  for  values  in  tables 

SHIFT 

Shift  index  values  (calculated  points  become  initial  points) 

VECTR 

Stores  calculated  values  in  one  array 

DNORM 

Calculates  normal  depth 

DCRIT 

Calculates  critical  depth 

WAV PLOT 

Main  plot  program 

SORTX 

Arrangement  of  array  for  plotting 

As  written,  the  program  has  several  limitations  which  may  require 
changing  for  different  problems.  The  only  end  boundary  conditions  pro¬ 
grammed  are  those  of  known  flow  rates  and  water  levels.  There  are  no 
provisions  made  for  varying  slope  or  cross-section  of  the  channel,  and  it 
is  also  assumed  that  the  roughness  is  constant  throughout,  regardless  of 
the  depth  of  flow.  Since  the  prediction-corrector  method  is  used  to 
obtain  the  solution,  accuracy  tests  are  required  and  these  may  need  chang¬ 
ing,  depending  on  the  accuracy  sought. 
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5 . 4  Testing  the  Program 

As  a  check  on  the  accuracy  of  the  program,  the  computer  solutions 
were  compared  to  the  field  observations.  Tests  were  initially  run  on 
open  channel  condi tions  and  then  the  program  was  tested  wi th  an  ice 
cover  imposed  on  the  channel.  From  the  available  field  data  the  channel 
cross-section  was  assumed  to  be  trapezoidal  with  a  base  of  210  feet  and 
side  slopes  1:15.  The  bed  slope  was  considered  to  be  the  average  slope 
over  the  portion  of  channel  analyzed,  while  the  Manning's  n  was  taken 
as  0.029  for  open  channel  flow  and  0.023  for  flow  under  the  ice  cover. 

For  the  open  channel  flow  tests,  a  section  of  the  North  Saskatchewan 
River  from  Lodgepole  to  Edmonton  was  used  as  the  test  section.  The  up¬ 
stream  boundary  condition  was  stage  versus  time,  while  the  downstream 
boundary  condition  was  the  rating  curve.  Surges  of  various  heights  and 
frequency  were  imposed  at  the  upstream  end  of  the  channel  and  they  were 
allowed  to  propagate  downstream  under  various  conditions.  Plots  of  velo¬ 
cities  and  depths  for  all  runs  performed  are  shown  in  Chapter  VI. 

The  tests  for  ice  cover  conditions  were  carried  out  with  some 
variations,  such  as,  the  test  section  used  was  that  of  the  North  Saskat¬ 
chewan  River  from  the  Brazeau  River  junction  to  Edmonton.  The  upstream 
boundary  condition  was  discharge  versus  time,  with  the  downstream  boundary 
condition  being  the  open  water  stage-discharge  curve.  The  imposed  dis¬ 
charge  at  the  upstream  end  was  the  sum  of  the  discharges  at  Rocky  Mountain 
House  and  the  Big  Bend  Plant.  The  fluctuations  in  flow  only  occur  at  the 
Big  Bend  Plant  gauge,  therefore  it  was  necessary  to  estimate  the  time 
required  for  these  variations  to  propagate  the  ten  miles  downstream  to 
the  junction  of  the  North  Saskatchewan  River.  The  remainder  of  the  test¬ 
ing  procedure  was  similar  to  that  of  the  open  channel  flow  condition. 


' 


CHAPTER  VI 


RESULTS 

6 . 1  Genera  I 

The  results  of  the  computer  simulation  of  wave  propagation  are 
best  shown  by  means  of  plots.  Plotting  depth  and  velocity  against  distance 
and  time  fully  describes  the  changes  which  occur  in  the  surge  as  it  pro¬ 
gresses  downstream. 

6.2  Plots 

The  plots  on  the  following  pages  are  the  results  obtained  from  the 
runs  performed.  As  the  channel  properties  were  varied  for  various  runs, 
the  conditions  imposed  on  each  run  are  listed  with  the  resulting  plots. 
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FIGURE  6-2  COMPARISON  OF  WATER  LEVELS  AT  DOWNSTREAM  END 
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FIGURE  6-4  COMPARISON  OF  WATER  LEVELS  AT  DOWNSTREAM  END 
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FIGURE  6-6  VELOCITY  AND  DEPTH  AT  DOWNSTREAM  END 
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FIGURE  6-8  VELOCITY  AND  DEPTH  AT  20.0  HOURS 
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FIGURE  6-10  VELOCITY  AND  DEPTH  AT  60.0  HOURS 
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FIGURE  6-12  COMPARISON  OF  Vv'ATER  LEVELS  AT  DOWNSTREAM  END 
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FIGURE  6-14  COMPARISON  OF.  WATER  LEVELS  AT  DOWNSTREAM  END 
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FIGURE  6-16  VELOCITY  AND  DEPTH  AT  DOWNSTREAM  END 
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FIGURE  6-18  VELOCITY  AND  DEPTH  AT  DOWNSTREAM  END 
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FIGURE  6-20  VELOCITY  AND  DEPTH  AT  20.0  HOURS 
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FIGURE  6-22  VELOCITY  AND  DEPTH  AT  60.0  HOURS 


CHAPTER  VII 


DISCUSSION  OF  RESULTS 


7. 1  Genera  1 

This  chapter  discusses  the  results  recorded  in  the  previous  chapter. 
The  graphs  provide  an  easy  means  of  comparing  the  various  results.  The 
discussion  mainly  deals  with  the  effect  of  channel  roughness  on  the  surge 
prof i 1 es . 

7 . 2  Open-Water  Surges 

The  results  of  open-water  surges  were  plotted  on  Figures  6-1  to 
6-10.  Both  the  calculated  and  observed  results  for  a  single  surge  were 
plotted  on  Figure  6-2.  The  surge  plotted  with  a  solid  line  represents  the 
calculated  results  for  the  conditions  most  closely  related  to  natural 
conditions.  The  travel  time  calculated  was  similar  to  the  observed  value 
but  the  observed  height  of  surge  was  approximately  25%  less  than  that 
calculated.  The  only  difference  between  the  two  calculated  runs  was  a 
7%  increase  in  the  value  of  n.  A  noticeable  flattening  of  the  surge  occurs 
as  a  result  of  the  increased  roughness.  The  decrease  in  surge  amplitude 
was  6%  while  the  increase  in  travel  time  was  6%. 

In  Figures  6-3  and  6-4,  an  increase  of  100%  in  side  slopes  and  a 
7%  decrease  in  n  resulted  in  an  increase  of  5%  in  travel  time  and  an  in¬ 
crease  of  20%  in  surge  height.  The  difference  in  initial  depths  at  the 
downstream  end  was  a  result  of  the  stage-discharge  boundary  condition. 

The  remaining  figures  of  open-water  surges  depict  a  series  of 
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surges  propagating  downstream  in  the  channel.  These  figures  show  predicted 
profiles  of  the  channel  at  various  times  as  well  as  the  predicted  stage 
recordings  for  various  stations.  This  gives  a  good  description  of  the 
surge  as  it  progresses  downstream.  The  computed  travel  time  was  8%  less 
than  the  observed  time.  The  amplitude  of  the  observed  surge  was  approxi¬ 
mately  15%  higher  than  that  calculated  for  the  largest  surge.  The  smaller 
the  amplitude  the  more  closely  related  were  the  calculated  and  observed 
values  for  both  travel  time  and  surge  height. 

In  general,  the  open-water  surges  behaved  as  expected.  The  varia¬ 
tions  between  predicted  and  observed  results  increased  with  an  increase 
in  surge  height.  An  increase  in  the  roughness  coefficient  increased  the 
travel  time  and  decreased  the  height  of  the  surge.  This  resulted  in  the 
surge  becoming  flatter  and  broader. 

7 . 3  Surges  Under  Ice  Cover 

The  remaining  figures  from  Figure  6-11  to  6-22  inclusive  are  the 
results  of  surges  propagating  down  a  channel  under  an  ice  cover .  In 
Figures  6-11  and  6-12  the  same  discharge  was  imposed  at  the  upstream  end 
but  the  conditions  of  the  channel  were  varied  for  each  run.  Figure  6  13 
and  6- 1 A  compare  surges  in  open  water  and  under  an  ice  cover.  The  remain 
i ng  runs  were  concerned  with  a  series  of  surges  propagating  down  the 

channel . 

An  increase  of  8%  in  the  value  of  n  produced  an  increase  of  1% 
in  travel  time  but  the  height  of  the  surge  was  reduced  by  only  a  slight 
amount.  The  shapes  of  the  surges  were  very  similar  in  both  cases.  An 
increase  of  10%  in  the  channel  slope  resulted  in  a  decrease  of  5^  In  the 
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travel  time  of  the  surge  and  an  increase  of  6%  in  the  height  of  the  surge. 

Therefore,  the  surge  had  a  sharper  crest  and  was  not  spread  out  over  as 
large  a  distance. 

Figures  6-15  to  6-22  are  the  results  of  tests  performed  on  a  series 
of  surges  propagating  down  the  channel.  The  only  difference  in  the  two 
runs  was  the  initial  discharge.  The  variation  in  initial  condition  of 
the  channel  only  effected  the  profile  of  the  first  surge.  The  channel 
stabilized  quickly,  resulting  in  subsequent  surges  having  the  same  profile 
regardless  of  the  initial  conditions  imposed.  The  observed  data  plotted 
on  the  figures  had  different  initial  conditions  than  the  calculated  results, 
thus  the  only  va 1  id  comparison  would  be  to  compare  the  second  surge.  The 
predicted  travel  time  was  12%  less  than  observed  time  and  the  predicted 
height  of  surge  was  36%  greater  than  observed  height  of  surge. 

From  the  tests  performed,  the  ice  cover  did  not  appear  to  have  as 
much  effect  as  observed  in  the  field.  There  was  a  large  variation  in  the 
calculated  and  observed  conditions  for  the  runs  under  ice  cover.  Although 
the  calculated  results  showed  that  the  surge  was  damped  considerably  due 
to  the  ice  cover,  it  was  not  as  large  as  the  actual  damping.  Possibly  the 
calculated  roughness  coefficient  was  not  as  large  as  the  actual  value. 

7 . k  Comparison  of  Open  Water  Surges  and  Surges  Under  an  Ice  Cover 

The  results  of  a  single  surge  propagating  down  a  channel  under 
various  conditions  are  tabulated  in  Table  7*1*  No.  2  and  are  the 
natural  conditions  for  the  North  Saskatchewan  River  and  are  compared 
in  Figures  6-13  and  6-1**.  For  open  water  the  decrease  in  surge  height 
over  the  reach  was  20%  and  under  an  ice  cover  it  was  **5$.  The  ice  caused 
a  considerable  damping  of  the  surge.  The  travel  time  of  the  surge  for 
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TABLE  7.1  -  SUMMARY  OF  COMPUTED  SURGE  DATA 


No. 

Channel 

Cover 

n 

S 

X10~3 

SS 

Ups  tream 
Surge 
Height 
(ft) 

Downstream 

Surge 

Height 

(ft) 

Travel 

Time 

(hrs) 

1 

Open 

0.027 

8.64 

30.0 

3.0 

2.8 

42 

2 

Open 

0.029 

8.64 

15.0 

3.0 

2.4 

41 

3 

Open 

0.031 

8.64 

15.0 

3.0 

2.2 

43 

4 

1  ce 

0.023 

8.64 

15.0 

3.9 

2.3 

41 

5 

1  ce 

0.025 

8.64 

15.0 

4.0 

2.2 

44 

6 

1  ce 

0.023 

9.50 

15.0 

3.9 

2.3 

39 

the  two  conditions  was  similar  even  though  the  initial  velocity  was  less 
for  the  ice  cover  condition.  This  occurs  because  the  surge  has  a  greater 
height  under  ice  cover  thus  the  celerity  of  the  surge  increases. 

Even  though  the  calculated  results  showed  some  large  fluctuations 
from  observed  values,  the  results  did  show  that  an  ice  cover  causes  con¬ 
siderable  change  to  the  profile  of  the  surge. 

7.5  Computer  Techniques 

The  computer  provided  a  very  necessary  tool  in  performing  the  cal¬ 
culations  for  the  numerical  analysis  of  this  study.  Discrepancies  in  the 
results  produced  by  these  techniques  are  small  compared  to  those  produced 
by  the  broad  assumptions  used  in  the  analysis.  In  Figures  6-15  and  o-17> 
the  discharge  data  was  very  irregular  but  small  errors  become  noticeable 
after  a  time  of  70  hours.  The  effect  produced  by  these  errors  was  negligible 
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CHAPTER  VIII 

CONCLUSIONS  AND  RECOMMENDATIONS 

8. 1  Genera  1 

A  review  of  the  literature  indicated  a  lack  of  information  on 
surges  as  well  as  the  roughness  coefficients  in  ice  covered  channels. 

The  present  study  was  initially  concerned  with  the  evaluation  of  composite 
roughness  coefficients  and  then  applying  these  results  to  surges  propa¬ 
gating  down  ice  covered  channels. 

8 . 2  Concl us  ions 

The  addition  of  an  ice  cover  results  in  considerable  change  in  the 
channel  roughness.  Even  though  the  numerical  coefficient  of  roughness  for 
an  ice  covered  channel  is  less  than  that  for  an  open  channel,  the  effect 
of  the  roughness  is  increased  as  a  result  of  the  increase  in  wetted  peri¬ 
meter.  During  the  mid-winter  period  on  the  North  Saskatchewan  River, 
Manning's  n  for  the  ice  cover  can  be  approximated  as  0.015-  This  is 
approximately  1/2  the  value  of  n  for  open  channel  conditions  at  low  flows. 
Thus  an  approximate  value  for  a  composite  roughness  coefficient  to  be 
applied  to  the  entire  channel  during  ice  cover  would  be  0.023.  On  the 
basis  of  the  roughness  coefficient  of  the  ice  cover  being  1/2  the  roughness 
of  the  bed,  the  depth  of  flow  would  increase  ]24  due  to  the  ice  cover. 

Stage  readings  would  increase  an  additional  amount  equal  to  the  bouyant 
thickness  of  the  ice  cover.  During  the  period  of  freeze-up  and  break-up, 
the  roughness  coefficient  is  increased  to  a  much  greater  extent  but  it  is 
very  difficult  to  estimate  this  increase  as  it  varies  from  year  to  year 
depending  on  the  climatic  conditions. 
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From  the  results  obtained,  the  ice  cover  had  a  significant  effect 
on  the  characteristics  of  a  surge.  The  ice  cover  had  little  effect  on 
the  travel  time  of  the  surge  even  though  the  flow  velocity  decreases. 

This  is  a  result  of  the  increase  in  celerity  of  the  surge  due  to  its 
increase  in  height.  The  ice  cover  also  causes  a  greater  damping  of  the 
su  rge . 

The  method  of  characteristics  combined  with  the  numerical  techniques 
applied  did  not  produce  as  accurate  results  as  would  be  desirable.  Much 
of  the  variations  in  results  can  be  contributed  to  the  limited  amount  of 
channel  geometry  applied  and  the  large  length  of  channel  considered.  As 
an  estimate  of  surge  properties  for  design  purposes,  this  method  could  be 
considered  since  it  gives  surges  at  the  downstream  end  which  are  higher 
than  observed  surges. 

8 . 3  Recommendat i ons 

In  view  of  the  present  study  there  is  a  great  need  for  more  study 
in  this  area  in  order  that  the  effect  of  an  ice  cover  on  surges  may  be 
better  understood.  It  is  therefore  recommended  that: 

(a)  there  is  further  study  on  the  composite  roughness 
coefficient  produced  by  an  ice  cover  with  particular 
attention  on  the  freeze-up  and  break-up  period. 

(b)  test  be  conducted  to  determine  the  effect  of  various 
changes  in  channel  geometry  and  boundary  conditions 
used  with  a  computer  program  in  an  attempt  to  verify 
or  limit  the  conclusions  of  this  study. 
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APPENDIX  A 


DATA  COLLECTED 
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TABLE  A. 1 


STAGE-DISCHARGE  READINGS  AT  EDMONTON 


DATE 

STAGE* 

DISCHARGE** 

DATE 

STAGE* 

DISCHARGE** 

OCT/69 

29 

2062.9 

2450 

1 

2064.1 

5190 

30 

2062.8 

2360 

2 

2063.8 

4360 

31 

2062.9 

2360 

3 

2063.7 

4080 

NOV/69 

4 

2063.7 

4160 

1 

2063.0 

2740 

5 

2063.7 

4270 

2 

2063.0 

2710 

6 

2063.7 

4270 

3 

2063.1 

2630 

7 

2063.5 

4010 

4 

2063.1 

2610 

8 

2063.4 

3800 

5 

2063.1 

2670 

9 

2063.4 

3740 

6 

2063.1 

2740 

10 

2064.2 

5440 

7 

2063.1  P 

2760 

11 

2063.6 

4080 

8 

2063.1 

2570 

12 

2063.6 

4020 

9 

2062.9 

2340 

13 

2063.6 

4080 

10 

2063.0 

2480 

14 

2063.6  F 

4060 

1 1 

2063.0 

2340 

15 

2063.4  P 

3850 

12 

2063.3 

2980 

16 

2063.4 

3580 

13 

2063.3 

2700 

17 

2063.4 

3480 

14 

2063.4 

2910 

18 

2063.3 

3340 

15 

2064.0 

4800 

19 

2063.1 

3180 

16 

2063.9 

4040 

20 

2063.0 

3010 

17 

2064.6 

3690 

21 

2063.0 

2980 

18 

2064.5 

1770 

22 

2063.3 

3450 

19 

2067.2 

5700 

23 

2063.1 

3110 

20 

2066.9 

4040 

24 

2063.2 

3240 

21 

2067.0 

2180 

25 

2063.1 

3120 

22 

2067.5 

3010 

26 

2063.1 

3050 

23 

2067.4 

3000 

27 

2063.0 

2780 

24 

2067.3  P 

3300 

28 

2062.9 

2610 

25 

2067.3  P 

3610 

< 
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TABLE  A. 1  CONTINUED 


DATE 

STAGE* 

DISCHARGE** 

26 

2067-2  P 

3150 

27 

2067.2 

3100  E 

28 

2067.3  P 

3100  E 

29 

2067.1 

3100  E 

30 

2066 . 6 

3100  E 

DEC/69 

1 

2066.4 

2670 

2 

2066.1  P 

2960 

3 

2066.0 

2890 

4 

2065.7  P 

2530 

5 

2065.8 

2700 

6 

2065.6 

2570 

7 

2065.4 

2470 

8 

2065.4 

2570 

9 

2065.1 

2340 

10 

2064.7 

2110 

11 

2064.7 

1980 

12 

2064.6 

1980 

13 

2064.5 

1990 

1  4 

2064.7 

2140 

15 

2064.9 

2470 

16 

2065.3 

2930 

17 

2065.1 

2720 

18 

2064.8 

2510 

19 

2064.1 

2330 

20 

2064.8 

2470 

21 

2064.7 

2520 

22 

2064.5 

2330 

23 

2064.4 

2320 

DATE 

STAGE* 

DISCHARGE** 

24 

2064.4 

2230 

25 

2064.4 

2260 

26 

2064.6 

2400 

27 

2064.7 

2620 

28 

2064.5 

2460 

29 

2064.5 

2440 

30 

2064.4 

2760 

31 

2064.4 

2760  E 

JAN/70 

1 

2064.3 

2480 

2 

2064.6 

2160 

3 

2064.6 

2410 

4 

2064.7 

2500 

5 

2065.0 

2790 

6 

2065.0 

3100 

7 

2064.9 

2810 

8 

2064.6 

2520 

9 

2064.4 

2270 

10 

2064.5 

2350 

11 

2064.4 

2290 

12 

2064.2 

2100 

13 

2064.2 

2050 

14 

2064.4 

2180 

15 

2064.5 

2240 

16 

2064.6 

2510 

17 

2064.6 

2620 

18 

2064.5 

2270 

19 

2064.4 

2220 

20 

2064.3 

2120 

' 

TABLE  A. 1  CONTINUED 


DATE 

STAGE* 

DISCHARGE** 

18 

2064.7 

2560 

19 

206*4 .3 

2830 

20 

2064 . 3 

2810 

21 

2064.3 

2790 

22 

2065.0 

2970 

23 

206*4.3 

2910 

24 

206*4.7 

2750 

25 

206*4.7 

2610 

26 

2065.0 

2710 

27 

206*4.3 

2700  E 

28 

206*4.3 

2700  E 

MARCH/70 

1 

2064.3 

2850 

2 

2064.3 

2930 

3 

2064.8 

2680 

k 

2064.7 

2570 

5 

2065.1 

2740 

6 

2065.0 

2590 

7 

2065.1 

2710 

8 

2065.5 

3340 

9 

2065.2 

2960 

10 

2064.3 

2570 

11 

2065.0 

2520 

12 

2065.1 

2760 

13 

2065.1 

2830 

14 

2065.2 

2840 

15 

2065.2 

2830 

16 

2065.2 

2790 

17 

2065.1 

2800 

DATE 

STAGE*. 

DISCHARGE** 

21 

2064.4 

2230 

22 

2064.5 

2400 

23 

2064.5 

2390 

24 

2064.6 

2420 

25 

2064.5 

2390 

26 

2064.4 

2380 

27 

2064.3 

2220 

28 

2064.3 

2260 

29 

2064.3 

2210 

30 

2064.3 

2210  E 

31 

2064.3 

2210  E 

FEB/70 

1 

2064.6 

2450 

2 

2064.7 

2460 

3 

2064.6 

2420 

4 

2064.5 

2320 

5 

2064.8 

2580 

6 

2064.7 

2480 

7 

2064.6 

2420 

8 

2064.6 

2390 

9 

2064.4 

2260 

10 

2064.6 

2390 

11 

2064.7 

2420 

12 

2064.7 

2510 

13 

2064.7 

2590 

14 

2064.8 

2620 

15 

2064.9 

2790 

16 

2064.9 

2930 

17 

2064.7 

2630 

' 
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TABLE  A . 1  CONTINUED 


DATE 

STAGE* 

DISCHARGE** 

1  1 

2067.6 

9130 

12 

2067.2 

7800 

13 

2065.6 

6760 

14 

2064.7 

6200 

15 

2064.5 

6000 

16 

2064.3 

5150 

17 

2064.2 

4860 

18 

2064.1 

4720 

19 

2064.1 

4720 

20 

2063.9 

4560 

21 

2064.0 

4800 

22 

2064.3 

5100 

23 

2064.4 

5460 

2  4 

2064.3 

5230 

25 

2064.1 

4740 

26 

2064.1 

4700 

27 

2064.1 

4760 

28 

2064.0 

4380 

29 

2064.2 

6260 

30 

2064.0 

6000  E 

DATE 

STAGE* 

DISCHARGE** 

18 

2065.1 

2540 

19 

2065.2 

2920 

20 

2065.1 

2680 

21 

2065.0 

2440 

22 

2064.9 

2320 

23 

2064.8 

2210 

24 

2064.8 

2220 

25 

2064.9 

2330 

26 

2064.9 

2270 

27 

2064.9 

2290 

28 

2065.0 

2450 

29 

2064.9 

2360 

30 

2064.9 

2200 

31 

2064.9 

2210 

APRIL/70 

1 

2065.0 

2470 

2 

2065.1 

2460 

3 

2065.5 

3040 

4 

2065.8 

3950 

5 

2066 . 1 

4500 

6 

2066.9 

5630 

7 

2067.8 

8210 

8 

2067.8 

9890 

9 

2067.7 

8700 

10 

2067.7 

9030 

*  STAGE  AT  MAYFAIR  PARK  GAUGE  (CITY  ELEVATION) 

**  DISCHARGE  AT  LOW  LEVEL  BRIDGE  (CFS) 

P  PUMPING 
E  ESTIMATED 
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TABLE  A. 2 
CURRENT-METERING 

ON  THE  NORTH  SASKATCHEWAN  RIVER  AT  MAYFAIR  PARK  GAUGE 

ON  FEBRUARY  7,  1970 


SECTION 

DEPTH 

VELOCITY 

SECTION 

DEPTH 

VELOCITY 

(fps) 

(fps) 

#1 

FROZEN 

TO 

#5 

4 1  3" 

1  .07 

D  =  100' 

BOTTOM 

D  =  175' 

4 1 0" 

1.37 

n 

T.D.  =  4,6“ 

3'9" 

1.47 

D  =  1 25 1 

I.T.  =  1.3' 

3 1 6" 

1  .56 

T.D.  =  3'3" 

NO 

W. 1 .  =  +0.5" 

3  1 3" 

1  .68 

1  .T.  =  1.85' 

VELOCITY 

3'0" 

1  .73 

W.  1 .  =-1.5" 

2'9" 

1.73 

2 1 6" 

1 .82 

n 

2 1 3" 

1 .82 

D  =  150' 

2  1  0" 

1.49 

T.D.  =  4'3" 

NO 

1.77 

I.T.  =  1.5' 

VELOCITY 

1  ’9" 

1.59 

W. 1 .  =  +0.5" 

1'5" 

1.29 

if'3" 

0.85 

#6 

5' 1" 

0.78 

D  =  162' 

3'9" 

1.19 

D  =  200' 

4'9" 

1.14 

T.D.  =  4'6" 

3 1 3" 

1.23 

T.D.  =5' 4" 

4'6" 

1.26 

I.T.  =  1.5' 

2 1 9" 

1.45 

I.T.  =  1.9' 

V3" 

1.37 

W. 1 .  =  +0.5" 

2'3" 

1.39 

W. 1 .  =0 

4 1 0“ 

I.36 

2*0" 

1  .32 

1 .36 

1 1 7" 

1.13 

3'9" 

1.35 
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TABLE  A. 2  CONTINUED 


SECTION 

DEPTH 

VELOCITY 

(fps) 

SECTION 

DEPTH 

VELOCITY 

(fps) 

#13 

5 1 0" 

0.61 

D  =  600' 

3'3" 

0.92 

D  =  550 1 

V9" 

1.02 

T.D.  =  3'9" 

3'0" 

1.14 

T.D.  =  5 1 3" 

V6" 

1  .09 

I.T.  =  1.7' 

2'9" 

1 .33 

I.T.  =  1.5' 

4'3" 

1  .36 

W. 1 .  =  +0.5" 

2 1 6" 

1.43 

W. 1 .  =  +1.0" 

Vo" 

1.59 

2'3" 

1 .54 

3'9" 

1  .64 

2  1 0" 

1 .34 

3*6" 

1.77 

1  '9" 

0.98 

3'3" 

1.89 

#15 

3*4" 

0.47 

3 ' 0" 

2.05 

D  =  628' 

3 1 0" 

0.90 

2 1 9“ 

2.23 

T.D.  =  3 1 7" 

2'9" 

1  .04 

2 '6" 

2.36 

I.T.  =  17' 

2 ' 6" 

1  .06 

2 1 3n 

2.27 

W.  1  .  =  -1.0" 

2 1 3" 

1 .05 

2 1 0" 

2.26 

2'0" 

0.98 

1  ’9" 

2.05 

1  '9" 

0.71 

1 1 7" 

1 .30 

#16 

FROZEN 

TO 

#14 

3 '6" 

0.79 

D  =  650' 

BOTTOM 

L _ 

D  -  DISTANCE  FROM  WEST  BANK 
TD  -  TOTAL  DEPTH  OF  WATER 
IT  -  ICE  THICKNESS 

Wl  -  WATER  LEVEL  WITH  RESPECT  TO  THE 
TOP  OF  THE  ICE 
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LINEAR  REGRESSION 

The  values  of  velocity,  x,  and  the  logarithm  of  the  depth,  y,  were 
fitted  with  a  straight  line  by  the  method  of  least  squares.  The  most 
probable  position  of  such  a  line 

x  =  a  +  by  -  -  - - (B-l ) 

is  such  that  the  sum  of  the  squares  of  the  deviation  of  all  points  from 
the  line  is  a  minimum.  For  n  observations  this  occurs  when! 

2 

3  “  Ey  Ex  -  EyExy  -  -  -  -  (b-2) 

nEy2  -  (ly)2 

b  =  nExy  -  ExEy  _____  (3-3) 

nEy2  -  (Ey)2 


As  a  test  on  the  exactness  of  the  fit,  a  correlation  coefficient  was 

( 

computed, 


9 


C  =  nExy  -  ExEy 

'TnEx2  -  (Ex)2][nEy2  -  (E^) 2] 


-  (8-4) 


MANNING'S  n  FOR  TOTAL  CROSS-SECTION 


Values  of  n  were  obtained  at  several  sections  across  the  cross- 
section.  From  these  it  was  then  necessary  to  obtain  a  n  for  the  entire 
cross-section.  The  channel  cross-section  was  divided  intom  sections 
with  a  measurement  for  each  section.  The  sum  of  the  discharges  of  each 
section  must  equal  the  total  discharge. 

m 

=  E  dQ. 

i  =  l  ' 


a 


(B-5) 


. 


The  sections  are  wide  compared  to  the  depth  and  the  bed  slope  is  constant, 
giving: 


1.486  /d\2/3  s,/2  db 


m 

E 

i  =  1 


n 


bd 


5/3 


m 

E 

i  =1 


Mi 

n . 


5/3 


where: 

Q.  -  total  discharge 
dQ.  -  discharge  for  each  section 

n  -  Manning's  n  for  the  entire  cross-section 

n.  -  Manning's  n  for  each  section 

b  -  sum  of  the  widths  of  all  sections 
b.  -  width  of  each  section 

i 

d  -  total  cross-sectional  area  divided  by  b 
d.  -  depth  of  each  section. 


sl/2  d.b.  -  -  (B-6) 

-  (B-7) 


. 

TABLE  B.l 


AREA  AND  DISCHARGE  CALCULATIONS 


Section 

b. 

i 

ft . 

d. 

i 

ft . 

b.d. 

i  i 

ft2 

^Mean 

ft/sec 

dQ.  =  b.d.-  VM 

1  1  1  Mean 

cf  s 

3 

6.0 

0.0 

0.0 

0.0 

0.0 

4 

12.5 

3.00 

37-5 

1.22 

45.8 

5 

19.0 

3.20 

60 . 8 

1.54 

93.6 

6 

37.5 

3.40 

127.5 

1.32 

168.3 

7 

50.0 

3.50 

175-0 

1  .32 

231 .0 

8 

50.0 

4.45 

222.5 

1.20 

267.0 

9 

50.0 

2.85 

142.5 

1  .32 

188. 1 

10 

50.0 

4.35 

217-5 

1.17 

254.5 

1 1 

50.0 

2.55 

127-5 

1  .94 

247.4 

12 

50.0 

3.15 

157.5 

2.34 

368.6 

13 

50.0 

3.75 

187.5 

1.71 

320.6 

14 

39.0 

2.10 

81 .9 

1.13 

92.5 

15 

25.0 

1.95 

48.8 

0.85 

41.5 

16 

1 1  .0 

0.0 

0.0 

0.0 

0.0 

=  Eb.d./Eb.  = 

i  i  i 


AVERAGE  DEPTH 


3.17  FT. 


' 
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TABLE  B . 2 

REGRESSION  ANALYSIS 


SECTION 

BED  BOUNDARY 

ICE  BOUNDARY 

X  =  A 

+  BY 

CORRELATION 

COEFFICIENT 

X  =  A 

+  BY 

CORRELATION 

COEFFICIENT 

A 

B 

A 

B 

4 

1  . 243*1 

0.6549 

0.975 

1 .4329 

0.2814 

0.999 

5 

1.5536 

0.7301 

0.991 

1 .8159 

0.4915 

0.989 

6 

1  .2533 

0.6669 

0.950 

1 .4458 

0.2337 

0.983 

7 

1.3127 

0.4876 

0.984 

1.4621 

0.4095 

0.956 

8 

1.1609 

0.5115 

0.985 

1 .3912 

0.3963 

0.965 

9 

1 .3413 

1 .0894 

0.989 

1 .6238 

0.5121 

0.931 

10 

1.0595 

0.5838 

0.979 

1 .4722 

0.4647 

0.984 

11 

1 .9841 

1.3575 

Cs| 

00 

• 

0 

2.4197 

0.6506 

0.980 

12 

2.2986 

1.7512 

0.993 

3.2528 

1.19711 

1  .000 

13 

1 . 4555 

1 .6342 

0.978 

2.4520 

0.9587 

0.955 

14 

1.3176 

1 .0010 

0.977 

1.6748 

0.6503 

0.997 

15 

1  .0232 

0.8175 

0.950 

1.1214 

0 .3668 

0.984 
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TABLE  B.3 
COMPOSITE  n 


n 

n 

COMPOSITE  n 

SECTION 

FOR  BED 
BOUNDARY 

FOR  ICE 
BOUNDARY 

WEIGHTED 

AVERAGE 

PAVLOVSKIY 

BELOKAN  S 
SABANEER 

LARSEN 

A 

0.0290 

0.0098 

0.0194 

0.0217 

0.0206 

0.0184 

5 

0.0255 

0.0140 

0.0197 

0.0206 

0.0202 

0.0201 

6 

0.0292 

0.0081 

0.0186 

0.0215 

0.0201 

0.0158 

7 

0.0200 

0.0146 

0.0173 

0.0176 

0.0174 

0.0173 

8 

0.0236 

0.0149 

0.0193 

0.0198 

0.0195 

0.0193 

9 

0.0476 

0.0167 

0.0321 

0.0358 

0.0340 

0.0302 

10 

0.0291 

0.0167 

0.0229 

0.0238 

0.0233 

0.0220 

11 

0.0392 

0.0137 

0.0265 

0.0295 

0.0280 

0.0266 

12 

0.0415 

0.0203 

0.0309 

0.0328 

0.0318 

0.0315 

13 

0.0580 

0.0219 

0.0400 

0.0440 

0.0420 

0.0419 

14 

0.0458 

0.0219 

0.0338 

0.0360 

0.0349 

0.0349 

15 

0.0528 

0.0174 

0.0351 

0.0395 

0.0373 

0.0307 
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TABLE  B.4 

MANNING'S  n  FOR  ENTIRE  CROSS-SECTION 


SECTION 

b.d.5/3 

i  i 

WEIGHTED  AVE. 

PAVLOVSKIY 

BELOKAN  S  SAB. 

LARSEN 

b.dT/3 

1  1 

b.d.3/3 

1  t 

b.d.5/3 

1  1 

b.d.5/3 

1  1 

n. 

1 

n. 

! 

n. 

1 

n. 

1 

4 

78.0 

4,020.6 

3,594.5 

3,786.4 

4,239.1 

5 

132.0 

6,700.5 

6,407.8 

6,534.7 

6,567.2 

6 

288.3 

15,500.0 

13,409.3 

14,343.3 

18,246.8 

7 

403.4 

23,317.9 

22,920.5 

23,183.9 

23,317.9 

8 

602.0 

3M91.7 

30,404.0 

30,871 .8 

31,191.7 

9 

286.4 

8,922.1 

8,000.0 

8,423.5 

9,483.4 

10 

579.6 

25,310.0 

24,352.9 

24,875.5 

26,345.5 

1 1 

238.0 

8,981.1 

8,067.8 

8,500.0 

8,947.4 

12 

338.4 

10,951.5 

10,317.1 

10,641.5 

10,742.9 

13 

452.6 

11,315.0 

10,286.4 

10,776.2 

10,801.9 

14 

13^-3 

3,973.4 

3,730.6 

3,848.1 

3,848.1 

15 

76.1 

2,168.1 

1 ,926.6 

2,040.2 

2,478.8 

WEIGHTED  AVERAGE 

n 

=  0.0225 

PAVLOVSKIY 

n 

=  0.0238 

BELOKAN  S  SABANEER  - 

n 

=  0.0231 

LARSEN 

n 

=  0.0219 
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C  WAV  21 

C  USES  SUBROUTINES 

C  INCLUDES  THE  CONDITION  OF  ICE  COVER 
REAL  MN 

COMMON  XO( 300,2 ) , VOl 300 ,2),T(251,2),X(251,2),V(251,2),Y(251,2), 
1H0G ( 20  ,4 )  , VECT ( 1004 ) , 3 , SO, XR , MN, TSCL , XSCL , SS, SI ,NNX< 2) 

1  WRITE (6, 1000) 

1000  FORMAT  Cl*  ,COX  ,  "WAV21  •  ) 

TSCL=1 ./3600. 

X  S  C  L  =  1  ./5280. 

R  EAD ( 5 , ICO  1 , END  =  999 ) B , XR , SO , MN , $$ 

1001  FORMAT! 5F15.6) 

S  I=SS  +  SQRT ( l.+SS*SS) 

WRITE(6,1002)B,XR,SO,MN,SS,SI 

1002  FORMAT! ' 06= • , F 10. 3, 5X , ' XR= ' , F10. 2, 5X , • S0=' ,F 10 .6 , 5X , ' MN= ' ,F10.3,5X 
1 , 'SS=' ,F10.2 , '  Sl=  ' , F 10 . 3 ) 

C  READ  IN  TABULATED  DATA 

C  K= 1  FOR  TIME  IN  HRS  ( XO )  VS  DISCHARGE  (YO)  AT  LEFT  END 
C  EACH  TABLE  PROCEEDED  BY  HOG  CARD  £  TERMINATED  WITH  X0=-9999.9 
K=  1 

7  NMX ! K ) =0 

READ! 5,  1003) (HOG! J,K)  ,  J  =  l,20) 

1003  FORMAT ( 20A4 ) 

DO  10  1=1, 3C0 

READ(5,1004  )X0(  I  ,K)  ,Y0(  I  ,.<) 

1004  FORMAT! 2F10. 3) 

IF( XO ( I , K ) .EQ. -9999. 9) GO  TO  12 
NMX!  K)=NMX( K)+l 
10  CONTINUE 
12  NX= { NMX (  K )  + 1 1 ) / 12 
WRITE! 6, 1015 )K 
1015  FORMAT ( 'OT ABLE ', 15) 

WRITE (6, 1003) ( HDG ( J , K )  ,J=1,20) 

I  2  =  0 

DO  20  1=1, NX 
11=12+1 
12=  I  1  +  11 

IF!  I2.GT.NMX(K)  )  I2  =  NMX(K) 

WRITE! 6, 10  05) (XO! J ,K) , J= 1 1 , 12) 

100  5  FORMAT  (  '0X0*.  '  ,  12F10.2) 

WRITE (6, 1006) ( Y0( J,K) , J=I1,I2) 

1006  FORMAT!'  YO : ' , 1 2F 1 0. 2 ) 

20  CONTINUE 

IF! K. EQ.2 ) GO  TO  23 

C  TRANSFORM  TIME  IN  HOURS  TO  TIME  IN  SECONDS 
K=NMX( 1 ) 

DO  22  1=1, K 

XO!  I  , 1 )  =  X0 (1,1 ) *3oOO. 

22  CONTINUE 

C 

C  READ  IN  STAGE-DISCHARGE  CURVE  AT  DOWNSTREAM  (RIGHT  HAND)  END 
C  K=2  FOR  DEPTH  IN  FEET  ( XO )  AND  FLOW  IN  CFS  (YO) 

C  CARDS  RRECEEDED  BY  HDG  CARD  £  TERMINATED  WITH  X0=-9999.9 
K=2 

GO  TO  7 

C  READ  IN  NO.  OF  INITIAL  POINTS  ON  REACH  £  CALCULATE  INITIAL  CONDITION 

23  READ! 5,  1007JN,  TMAX 
100  7  FORMAT ( I  10, F 10. 2) 

WRI TE ( 6, 1013 )N, TMAX 

1013  FORMAT (  ' 0N= '  , I  5 , 5X  ,  '  T  M  A  X  = ' , F20. 3,  '  HOURS') 

C  STORE  B » X  R  >  MN , SO , N , S  S ,  £  TMAX  IN  FIRST  TAPE  RECORD  SPACE 
VECT!  1  )=B 
VECT ( 2  )  =  XR 
VEC  T ( 3 ) =S0 


■ 

.  . 
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V  ECT ( 4 )=MN 
VECUS)  =  SS 
VECT(6)=10*N 
VECT  (  7  )  =  Tf'  AX 
WRITE (3)  VECT 
XR  =  t>2  80.*XR 
TMAX=TMAX*36C0 . 

I  F  (  N .  G  T  •  25  0  )  N= 2  5  0 
DX=XR /FLOAT ( N- 1 ) 

C  INITIAL  CONDITION  IS  THE  NORMAL  FLOW  E  DEPTH  CORRESPONDING  TO  THE 
C  INITIAL  UPSTREAM  DISCHARGE 

QNGRM= YO  ( 1 »  1 ) 

CALL  DNORN! QNORM ,  KN , B, SO* SS » YNURM ) 

CALL  DCR I T ( QNORM , 3  ,  SS , YCR I T) 

W  R 1  T  £  (  6 , 1 0  0  8  ) YNORM , YCR I T , QNORM 

10U8  FORMAT ( ‘ ONGRMAL  DE  PTH= •  ,F10.2,10X, 'CRITICAL  DEPT H= * , F 1 0 . 2 , 5X , • AT  Q 
I  =  '  ,  F 10 • 2  ) 

C 

C  SET  UP  INITIAL  CONDITIONS 
ARE  A=  YNORM"-  (  B+YNGRM*SS  ) 

Vl= QNORM/ (AREA) 

Yl= YNORM 
DO  30  1  =  1,  N 
T( I , 1 ) =0.0 

X  (  1,1) =DX*FLOAT( 1-1) 

VI  I , 1 ) = V 1 
Y(I,l)=Yl/ 

30  CONTINUE 

TlN+1, 1)=-1.0 
CALL  VECTR(N) 

WRITE! 3)  VECT 
I PNT  R= 1 

WRITE  (6,1009) 

1009  FORMA  I ( ' 0  INITIAL  CONDITIONS*) 

CALL  SCRIB 
WRITE (6, 101 4) 

1014  FORMAT! '  CALCULATED  RESULTS*) 

CALL  SORT  ( 1 PNTR, YKAX, VMAX, VMIN) 

C 

33  CALL  REACH! 1, IT) 

CALL  SHIFT 

C  VECTR(N)  TRANSFERS  T,  X,  V,  E  Y  TO  VECT  E  CHANGES  T  FROM  HOURS  TO 
C  SECONDS  AND  X  FROM  FEET  TO  MILES. 

CALL  VECTR(M) 

W R I T E  (  3  )  VECT  .. 

I  PNT  R=  I  PNT  R  i-l 

CALL  SORT  !I PNTR, YMAX, VMAX, VMIN) 

TM=T (1,1) 

XM=X ( 1,1) 

VM= V (  1,1) 

YM= Y { 1,1) 

CALL  BDY 1  (  T M ,  X M ,  V K ,  Y .'1 ,  T N  ,  X N  ,  V M ,  Y N ) 

T (  1 ,2)  =  TN 
X (  1 , 2  )  =  XN 
V!  1,2) =VN 

Y  (  1 »  2 )  =  Y  N 

CALL  REACH! 2 , IT) 

TL=T{  I  T, 1) 

X  L  =  X  {  I  T  ,  1 ) 

VL  =  V (  I  T  ,  1  ) 

YL= Y ( I T, 1) 

CALL  BDY2  ( TL , XL , VL  , YL , TN , XN , VN , YN ) 

I  T  =  I T  + 1 
T (  IT , 2  )  =  TN 
X! IT, 2 )=XN 
V ( IT , 2 ) =  V N 


■ 


. 
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Y(  IT, 2)  =  YN  u  n 

T( IT+1,2)=-1.0 
CALL  SHIFT 
CALL  VECTR(N) 

WRITE (3)  VECT 
IPNTR=IPNTR*1 

CALL  SORT  (  I  PNTR ,  YMAX  ,  VMAX ,  VM  IN  ) 

C  FIND  SMALLEST  T  CALCULATED  IN  LAST  PASS 
TST= . 996  1U 
DU  40  1=1,251 
I  F { T (  1,1)  .LT.O.  ) GO  TO  41 
I  F { T {  I  ,  1  )  . L r. TST) TST=T( 1,1) 

CONTINUE 

* 

IF ( TST  .  LT . TMAX ) GO  TO  33 
* 

TX=TMAX*TSCL 
WRITEI6, 1010)TST,TX 

FORMAT ( • 0  TST=  ’  ,  F 10 • 2 , 5X  ,  '  TMAX  =  * , F 10 .2 , 1  -END  OF  XT  PLANE  CALCLN '  ) 
WRITEI6, 1020) 

FORMAT ('0  MAXIMUM  AND  MINIMUM  VALUES  FOR  STEP  l1,/) 

WR I TE  < 6, 1021 )  Y M A X ,  VMAX,  VMIN 

FURMAT ( *0' ,  10X ,  *  YM AX :  '  , Fo. 2, 10X,  ' VMAX: ' , F 8 . 2 , 1  OX , ' VM 1 N : ' , F8 . 2/ ) 

C  FILL  LAST  RECORD  WITH  -1.0  VALUES 
DO  51  1=1, 1004 
VECT { I ) =  —  1 . 0 
51  CONTINUE 

C  OVERWRITE  VECT (1000)  TO  VECT (1002)  WITH  YMAX ,  VMAX,  VMIN 

VECT (  1000 )  =  YMAX 
VEC  T ( 1001 )  =  VMAX 
VECT { 1002 )=VM I N 
I PNTk= IPNTR+1 
WRITE (3)  VECT 
999  WRITE(6,1012) IPNTR  - 

1012  FORMAT (•  STEP1  COMPLETE • , I  20 , •  RECORDS  ON  TAPE') 

CALL  EXIT 

STOP 

END 

SUbROUT  I  NE  8D  Y  1  (  TM ,  XM ,  VM  ,  YM  ,  Tfl ,  XN  ,  VN  ,  YN) 

C  FUR  TABULATED  FLOW  AT  LEFT  HAND  END  ( N=  1  )  WITH  ICE  COVER 
C  TRAPEZOIDAL  CHANNEL  WITH  CONSTANT  SLOPE 

C  REQUIRES  SUBRIUTINE  TA3L4  WHICH  RETURNS  TWO  VALUES  OF  Q  UN  EITHER 
C  SIDE  OF  ESTIMATED  TN. 

C  THIS  SUBROUTINE  CALLED  BOY  I LQ  IN  NOTES 
REAL  MN 

COMMON  X0( 300,  2 ) , Y0(300,2 ) ,T( 25 1, 2 ) ,X{ 251 ,2) , V( 251 ,2) , Y( 251 ,2)  , 
1HDG ( 20 , 4 ) , V  E  C  T (  1004)  , 3 , SO , XR , MN,TSCL , XSCL , SS , S I ,NMX(2) 

C 

C  INITIALIZE 
B0  =  B 
XN=0. 

YN=YM 
VN  =  VM 

AM  =  Yi*1*  ( BO  +  YM-S  S ) 

CO M= SORT (32 . 1 7* AM/ ( B0+2 • -YM-SS) ) 

TN=TM-XM/ (VM-COM) 

5  F  M  =  V  M  *  V  M  -  M  N  *  M  N  -  (  2  .  *  (  DO  *-  Yr S  l  )  /  (  YM-  (  BO  +  YM-SS  )  )  )  -*  1  .333/2.22 
CALL  TABL4  (  1 , TN ,TM2 , Q  M2 , TM1 , QM 1 ,T  PI , OP  1 , TP  2, QP2 ) 

C  CURVES  SMOOTHED  BY  TAKING  AVERAGED  POSITIONS 
QM2= (QM1+0M2 ) -0.5 
T M2  =  (  TM1KM2  )  *0  . 5 
0  M  1  =  (  Q  M 1  +  Q  P 1 )  -  0 . 5 
T  M 1  =  (  T  M 1  +  T  P  1  )  -  0 . 5 
QP1=( QP1*QP2) *0.5 
TP  1 = ( TP1+TP2J-0.5 
QP2=QP2+QP2-QP1 


40 

C 

41 

C 


1010 

1020 

1021 


■ 

- 
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TP2=TP2+TP2-TP1  C_5 

ALF A=QM2 / (  (  TM2-TiU  )  *  (  TM2-TP1 )  *  {  TM2-TP2  )  ) 

BET  A  =  Q  1 1  /  {  (  TM1-TM2  )  *(  TM1-TP1 )  *(  TM  1-TP2  )  ) 

GAHA-  0  Pl/{  (  TP  1-  TM2  )  *  (  TP  l-T.'ll }  *  (  T  P1-TP2  )  ) 

DELT  =  QP2/(  (  TP2-TM2  )*(  TP2-TM1  )*(  TP2-TP1  )  ) 

AO=-ALFA*  TM1*T  P  1  ^ T P2i  —[it  I  A*TM2*TP1*TP2-GAMA*TM2  *TM  1*  T  P2- Dl  LT  *TM2* 
l  TM1*TP1 

A 1 = ( GAMA+DEL f ) *TM2*TM1  +  ( BETA  +  DEL T ) * T M2* TP  1  +  ( BE TA+GAMA ) * TM2* TP2 
1  +  (ALF  A  +  DEL  T  )  *T  Ml*  TP1  +  (  ALF A+GAMA  )  *  TM  1  *  ]  P2  +  (  ALF  A*-  3  E  T  A  )  *  T  P  1*  TP  2 
A 2= !-BETA-GAMA-DELT)*TM2+(-ALFA-GAMA-DELT)*TMl+ ( - ALF A-BET A-OELT )  * 
1TP1 +( -  ALFA-BET  A -GAM A)* TP  2 
A3= ALF A  +  BE  TA  +  GAMA +DELT 

C  NOTE  THE  EQUATION  Q.\=  AC  +A L  T N +A2*TN*TN+  A 3*TN*TN* TN  USES  THE  SAME 
C  COEFFICIENTS  FOR  ALL  ITERATIONS 
I TRN=0 
TTS  T  =  1  . 

VTST=1 . E-04 
YTST=l.E-04 
C 

10  QN=( ( A3*TN+A2 ) *TN+A1 )*TN+AO 
IFIQN.LT.O. ) WRITE (6, 1002) TN,QN 
1002  FORMAT ( *  BDYI,  TN= • , E15 .6 , •  QN=',E15.6) 

AN=YN*( B0+YN*SS) 

VN=QN/ AN 
PN=2.*( BO+YN*S I ) 

BN=B0+2.*YN*SS 

SFN  =  1 . 13506*QN*QN*HN*MN* I BO+YN*S l ) **1 . 333/  {  ( YN* ( BC+YN*SS )  ) **3.333) 
C ON  =  SORT ( 32 . 1 7* AN /BN) 

I TRN= I TRN+ 1 

C 

F 10=0 . 5* ( COM+CON ) * ( VN-VM ) -32 . 1 7* ( YN- YM ) +8 .0425* ( COM+CON ) * ( SFM-2.* 

1 SO+SFN )* ( TN-TM ) 


F  2  0 = 0 . 5  *  ! T  N  ■ -  T  M )  *  ( V  M  -  C  0  M  +  V  N  • -  C  0  N  )  +  X  H 

OF  1 0 TN  =  0 . 5  * ( C  OM+C  ON ) * ( 1 6 . Ob 5 * ( SFM-2 . *S 0  +  S FM )  +  { 1 . /AN+3  2. 17 * ( TN-TM ) * 
1SFN/QN)*<A1+2.*A2*TN+3.*A3*TN*TN) ) 

OF  ID YN =-3  2.  17+4.02125*1  1.-2 .*SS*AN/ ( BN* BN)  ) *-( 2. *( VN-VM) +32. 17*( TN- 
1 TM ) * ( SFM-2 . *  S0  +  SFN )  ) /CON  +  0 . 65 1* ( TN-TM)*! C OM  +  C ON ) *QN*Q N*  MN *  MN *  PN  *  *  0 


1 .333*1  SI-1 .2  5*Pi\*dN/AN)  / { AN**3 . 333 )-QN*BN/ l  2  .  *Ai\*  AN  )  *  l  COM+CON ) 
DF  2DTN=0 • 5* { VM-COM  + V  N-CON+ ( TN-TM)*! A1+2.*A2*TN+3.*A3*TN*TN ) /AN) 
DF2DYN=0. 5*( TM-TN ) * { VN* BN/ AN+ 1 o . 0 85* ( 1 .-2 .*SS*AN/ ( BN* BN) ) /CON) 
DEN=DF 1DTN*DF2DYN-DF 1 DYN*0F2DTN 
DT= { F 2 0*DF 10YN-F 1 0*DF20YN) /OEN 
DY= (F10*DF2DTN-F20*0F1DTN)/DEN 
T  N  =  T  N  +  0 . 5  *  D  T 
YN= YN+O . 5*0Y 
I F  {  ITRN.GT . 40 ) GO  TO  17 

1 F { (DT.GT.TTST ) .OR. (DY.GT . YT ST ) ) GO  TO  10 
GO  TO  20 

17  WRITE(6#100l)l TRN 

1001  FGRMATCO*  ,110,  '  ITERATIONS  IN  BDYI,  STOP') 

STOP 

20  QN= ( ( A3* TN+A2 ) *TN  +  A1 )*TN+AO 
V  N =  Q N / ( YN*BO  +  YN*YN*SS) 

TTN=TN/3600. 

W  R I T  E ( 6 , 100  3) T  TN , QN 

1003  FORMAT! *  UPSTREAM  END,  T=',E15.6,'  Q=',E15.o,'  C.F.S.') 

RETURN 


C 

C 

C 


END 

SUBROUTINE  REACH  (I  I, IT) 

CALCULATES  INTERIOR  POINTS  ON  X-T  DIAGRAM 
II  IS  INDEX  NUMBER  OF  FIRST  POINT  CALCULATED 
IT  IS  INDEX  NUMBER  OF  LAST  POINT  CALCULATED 
REAL  MN 

COMMON  X  C(  3  00, 2  )  ,  Y0(  300 ,2  )  ,  T  (  25  1 , 2  )  ,  X  <  251 ,  ) 

1HDG ( 20 , 4 ) , VECT ( 1004 ) , B , SO , XR , MN , TSCL , XSCL , SS 

K= I  I - 1 


V ( 25 i ,2) ,Y<  451,2) , 

SI ,NMX( 2) 


■ 

* 

' 

* 


C-6 


G  =  32. 17 
DO  100  1=1,250 
M  l  T  =  0 

1F( T( I +1 » 1 ) .LT.-O. 1 ) GU  TO  110 
TL=T< 1,1) 

XL -X { 1,1) 

VL  =  V( 1,1) 

YL  =  Y  (  1,1) 

TM=T( 1+1  ,1) 

XM=X( 1+1,1) 

VM=V( 1+1,1) 

YM=Y( 1+1,1) 

AL  =  YL  * ( B  +YL  *SS ) 

PL=B+2 .  *YL*SORT ( 1 . +SS*SS ) +B+2.*YL*SS 
SFL=  VL*VL*H.\I*MN/  (  2 . 22-  (  AL/PL)  **  1  .3333  ) 

AiY=  YM*  (  B+YM-SS  ) 

PM=B+2.*YM*SQRT( l.+SS*SS)+B+2.*YM*SS 

SFM=  VM*ViM*HN*MN/ (2. 22* ( AM/ PM) ** 1 .3333) 

C 

X  N  =  0 . 5  *  (  X  L  +  X  M  ) 

VN=0.5*( VL+VM) 

YN-0 .  5*( YL+YM) 

CLO=SQRT (32.17*AL/ <  S+2 .*YL*SS) ) 

CMO=SQRT (32. 17*AM/{ 3+2.*YM*SS) ) 

TN=TL+ ( XN-XL ) / ( VL+CLO) 

C 

10  Ai\|= YN* ( B+YN*S  S ) 

PN=2.*B+2.*YN*SQRT  (  1  .+SS*SS  )  +  2.*Yl\l*SS 
BN=B+2 .*YN*SS 

SFN-VN*VN*-1N*MN/2 . 2  2*  (  PM /AfJ ) ** 1 . 3333 
CNG=$QRT  (  G  *  AN/  (  B+2  .*YN*SS)  ) 

CUP=CLU+CN0 
COM=  CMG+CNO 
ISW=i 

F1T=0.5*(VL  +VN+COP ) 

F1X=-1.0 
F1V=0.5*( TN-TL) 

F 1 Y=G/ ( 4 .  *CNO ) * ( TN-T  L ) * { 1.-2.* SS* AN/ ( BN*BN) ) 

F2T  =  0.5*(VM+VN-CCJM) 

F2X=-1.0 
F2V=0.5*(TN-TM) 

F2  Y  =  -G/ ( 4 . *CNO ) * ( T  N-TM ) * ( 1 . -2 . * S S* AN/ ( BN*BN ) ) 

F3T=G/4. *C0P* ( SFL+SFN-2 .*$0 ) 

F  3V  =0 . 5  *  C  0  P  +  G  *  C  U  P  /  4 . 4  4  *  (  T 1 1  -  F  L  )  *VN*MN*MN*  (  PM/ AN)  **1.33  33 
F3Y=(0.5*(VN-VL )  +  3/4.* (TN-TL ) *( 5FL-2 .*SO  +  SFM)  ) *  C  G / ( 2. *CNO)  ) *(  1.-2. 
l*$S*AN/(  BN* BN )  ) +G+G/3. 33 *( TN-TL ) *CUP*  VM* VN*MM*MN* ( Pl\/AN)**0 . 3333* ( 
2  S I /AN- PN*8N/ ( 2 • *AN*AN) ) 

F4T=G/4. *C0M* ( SFM-2 . * S0+ SFN ) 

F 4  V  =  0 . 5*C0M  +  G/4.44*CCiT*  (  Tim-TM  )  *  VN*MN*  MN*  {  PN/ AM  )  **  1 . 33  33 

F4Y  =  (  0 . 5  *  (  VM-Vn  )  +G/4.*  (  TIM- TM  )  *  (  S  FM-2  .  *  SO  +  SFN  )  )  *  (  G/  (  2  .  *LNG  )  )  *  (  1 . -2  . 

1  *  S  S  *  A  N / ( BN* BN )  ) -G+G/3 . 33* ( TN-TM) *C0M*VN*VN*MN*MN* { PN/AN )**0 . 333  3*  ( 

2  S I / AN-PM*BN / { 2 . *aN*AN ) ) 

F  10  =  - 1  .  * ( 0 . 5  * ( TN-TL) *( VL+VN  +  COP) -(XN-XL) ) 

F  20=-  1  .  * ( 0 . 5  * (  TM-TM) *( VM+VN-COM) - ( XN-XM) ) 

F 3 G  =  -  1  . * ( 0 . 5 *C 0  P* ( VN-VL ) +G* ( YN-YL )  +G*0. 5*C0P* (  ( SFL  +SFN ) *0 • 5- SO) *1 T 

1N-TL) ) 

F40  =  -l.*(0. 5*C0M*  (  V.N-VM  )  -G*  (  YN-YM)  +G*0. 5*C0M*  (  (  SFM  +  SFM)  *0.5-S0)  *(  T 
1 M-T  M ) ) 

F 1 1 0=  F2  0-F  1  0 
FI  1T  =  F2T-F1T 
F 1 1 V  =  F2  V-F 1 V 
FI  1Y  =  F2Y-F  LY 

DET=F11T*(F3V*F4Y-F3Y*F4V)+F11V*(F3Y*F4T-F3T*F4Y)+F11Y*(F3T* F4 V-F 3 
1 V  *  F  4  T  ) 

DEL  T  =  F 1 10* ( F  3V* F4 Y-F 3Y*F4V )  +  F 1 1 V* ( F3Y*F40-F3  0*F4Y) + F 1 1 Y* ( F30* F4V-F 
1 3V  *F40 ) 


. 


' 


o  o 


C-7 

.  ^^7l:}l™F30^F4Y-F3Y^^)+F110«(F3Y^F4T-F3T^F4Y)+FllY«(r-iT^f40-F 

1  jU  ^  J-  4  T  ) 

1  ,^1;^;FJ1T:i:(F3V;!cF"t'J-F30*F4V)  f  FI  1V*(F30*F4T-F3T*F40)+F110*(F3T*F4V-F 

1  j  V  F  4  I  J 
DT  =  DEL T/OET 
DV=DELV/D£T 
DY=DELY/D£T 
T 1=TN+DT 
V1=VN+DV 
Y 1= YN+OY 

X 1 -XL  +  (  (Tl-TL  )/2.  )^(Vl  +  CLO  +  Vl  +  SQRT  (G-«iYlJMB+Yl*SS)/(P  +  2.*Yl*SS)  )  ) 
MIT=MIT+1 

IF  (MIT  .LT.  IQ)  GO  TO  20 
W  R I T  E  {  6  , 1 0  G  0 )  MIT, I 

-■FIT  E  (  6 , 100  1 )  TN  ,  XN  ,  YN  ,  VN  ,  T  1 ,  XI  ,  Yl ,  VI 

WR I TE ( 6 , 1001 ) FI  0 , F20, F30  , F40 » OET, DELT, QELV,DELY 

1000  FORMAT  (2X»  1  !■=•  ,  I  5 , 1  OX  ,  •  M I  T=  •  ,  I  5 ) 

1001  FORMAT ( IX, 8  E 1 4 . 6 ) 

20  GO  TO  (30,40,50,60) , ISW 
30  TST=ABS( Tl-TN)— 1.0 
GO  TO  70 

40  TST=ABS(X1-XN)-10.0 
GO  TO  70 

50  TST=AGS( Vl-VN) -1 .OE-4 
GO  TO  70 

60  TST=ABS( Y1-YN)-1 .OE-4 
70  IF( TST.GT.O. )G0  TO  80 
ISW=I$W+1 

IF (  ISW.LT.5)  GO  TO  20 
80  TN=  T 1 
XN=X1 
YN  =  Yl 
VN= V  l 

I F (  ISW.LT.5JG0  TO  10 
K  =  K  +  1 
T ( K , 2 ) =TN 
X{ K,2)=XN 
V ( K  ,  2 ) =VN 
Y ( K , 2 ) =YN 
100  CONTINUE 
C 

110  I T=K 

T { I T+ 1 , 2 )=- 1 . 0 

RETURN 

END 

SUBROUT  I NE  BDY2 { TL , XL , VL , YL , TN , XN , VN , YN) 

C  FOR  TABULATED  STAGE-DISCHARGE  CURVE  AT  kIGHT  HAND  END  OF  CHANNEL 
C  (DOWNSTREAM  END).  TRAPEZOIDAL  CROSS-SECTION  WITH  CONSTANT  SLOPE. 

C  REQUIRES  SUBROUTINE  TABL4  WITH  N  =  2. 

REAL  MN 

COMMON  XC(300,2),Y0(30G,2)  ,T(251,2)»X(251»2)  , V ( 2 5 1 , 2 )  ,Y(251,2), 

1HUG ( 20 , 4 ) , VECT ( 1004 ) , 3, SO, XR , MN , I SCL , XSCL , SS, SI , NMX( 2) 

INITIALIZE 
BO  =  B 
XN  =  XR 

AL=YL#(80+YL*SS) 

BL=B0+2.*YL*SS 
PL=2.*(30+YL*Si  ) 

COL=SQRT(32.17*AL/BL) 

SFL=VL*VL*MN*MN* ( PL/ AL ) ** 1 .33 3/2.22 
TN=TL+( XN-XL) / ( VL* COL) 

YN=  YL 

CALL  TABL4(2,YN,YV2,Q’12  ,YM1 , 0M1 ,  Y  PI  ,0P1  ,YP2,QP2) 

ALFA  =  QM2/(  ( YM2-YM1 ) *( YM2-YP1 )* ( YM2-YP2 ) ) 

BET A=QMl/( ( YM1-YM2 ) * ( YM1-YP1 )*( YM1-YP2 ) ) 


. 


,  -■  «  I  ifisfg  1  I  g  J- 


. 


GAMA=Q PI /< { YP 1- YM2 ) * ( YP 1-YM1 ) * { YP1-YP2 ) ) 

DELT  =  QP2  /  (  ( YP2-Y.M2  )  *  {  YP2-YM1 )  *  (  YP2-YP1 )  ) 

A0=-ALFA*YM1*Y  P1*YP2-BET  A*YM2*  YP  1*Y  P2-GAMA  *YM2  *YM  i  *Y P2- 0  EL  T*  YM 2* 
1YM1*YP1 

Al  =  (  GAMA+DELT  )  *  YM2*YM1+  (  BETA+DEL  T  )  * YM2*YP1  +  {  BE  TA+GAM A  )  *  YM2*YP2+ 

1 ( ALFA+OELT) *YM1*YP1+  ( ALF A+GAMA  )  *YM 1* Y P2  +  ( ALF A  +  b ET A ) * YP 1*Y P2 
A 2= ( -BETA-GAM A- DEL T ) *YM2  + ( -ALF A-GAMA-DELT )*YM1 + (-.aLFA-BET A-DELT ) * 
i YP 1  +  ( -AL  FA-BcT  A-GAMA ) *YP2 
A3=ALFA+BE  TA  +  GAMA  +  DEL  T 
C 

C  NOTE  THE  EXPRESSION  QN  =  AO  +  A  1*  YN  4-2  .  *YN**2  +A3  *YN**3  IS  USED  WITH  THE 
C  SAME  COEFFICIENTS  FOR  ALL  ITERATIONS. 

I  TRN  =  0 
TTS  T=1 • 

YTS  T= 1 .  E-04 

10  QN={ I A3*YN+A2 ) *YN+A1 )*YN+AO 
AN=YN*(BO+YN*SS) 

BN- 80  +  2 .  *YN*SS 
PN=2.*(B0+YN*S I ) 

SFN=  1 .  1350b*QN*QN*MN*MN*  (  B0+ YN*S  I )  **  1 . 333/  (  (YN-  (  80+YN*SS )  )  **3 .333 ) 
C0N=SQRT(32.17*AN/0N) 

VN=QN/AN 
I  TRN  =  I  TRN+ 1 
C 

0F1DTN=8 .0425* { CCL+CON) * ( SFL-2 . *S0 +SFN  ) 

DF2DTN=G.5*( VL+COL+VN+CON ) 

C 

OF  1 DYN  =  C . 5  *( COL  +  CON ) * ( AN* ( Ai +2 . *A 2*YN  +  3 . *A3* YN* YN ) -QN  *3N ) / ( AN*AN )  + 
132. 17+(C.5*( VN-VL )+b.0425* ( TN-TL )*< SFL-2 .*S0+SFN ) )*io .09*( 1 .-2.*SS 

2  * AN/ ( BN*3N )  )  / C  ON  +  J • 6 6 0* ( TN  — T L ) * ( C  OL  +  C ON ) *0 N* MN*  MN * PN*  * 0 • 3 33— ( S I +0 • 

3  75*PN* { 3 . * A3*YN*Y  N+2 • * A2*YN  + A1 ) -1 . 25*PN*QN*BN/AN) / ( AN**  3.33  3) 

C 

DF2DYN=C .5* ( TN-TL ) * (  ( A  1  +  2 . * A 2* YN +3 . * A3* YN* YN ) / AN-QN* BN/ (AN* AN)  + 
116.085*( l.-2.*SS*AN/(SN*BN) ) /CON) 

C 

F  10=C  .  5*  { COL  +  CON )  *(  VN-VL  )  +3  2. 17*(  YN-YL  )  +  8 . 04  25*  (  COL  +  C ON )  *  (  SFL-2  .* 

1  S 0  +  S  F  N  )  *  (  T  N -T  L  ) 

F2G=0.5* (TN-TL ) * ( VL+COL +VN+C ON) -XR+XL 
C 

DEN=DF 1 DTN*DF2  OYN-OF 1 0YN*DF2  DTN 
DT  = ( F2  0*DF1DYN-F 1 0*DF2DYN ) /OEN 
DY= ( F10*DF2DTN-F20*DFlDTN)/DtN 

C 

TN=TN+0. 5*DT 
YN- YN  +  0 . 5  *DY 
I  F  (  ITRN.GT  .  24  )  GCJ  TO  17 

IFKDT.GT.TTST)  .OR.(OY.GT.YTST)  )GU  TO  10 
GO  TO  20 

17  WRI TE ( 6, 1001 ) ITRN 

1001  FORMAT ( 1 0*  ,  1 10 » •  ITERATIONS  IN  BDY2*  *STOP**) 

STOP 

20  QN=( ( A3*YN+A2 ) *YN+A1 ) *YN+A0 
AN=YN*(B0+YN*SS) 

VN=QN/ AM 
RETURN 
END 

SUBROUTINE  TABL4 ( N , X P  ,  XM2 , YM 2 , XML  ,YM1 , XP 1 , Y P 1 , XP2 , YP 2 ) 

C  FINDS  TWO  PAIRS  OF  POINTS  ON  EITHER  SIDE  OF  TN 
C  TAKES  DATA  FROM  TABLE  N .  NM ( N )  IS  NUMBER  UF  TABLE  ENTRIES 
C  RETURNS  WITH  TWO  POINTS  ON  EITHER  SIDE  OF  XP .  EXCEPT  NEAR 
C  WHERE  THE  FOUR  END  POINTS  ARE  RETURNED  IF  TN  INSIDE  RANGE 
REAL  MN 

COT  MON  XC(300»2)»V0(3  00j2)»T(25l»2)jX(251»2)»V(251/2)»Y(25l»2)» 
1HDG(  20  »  4  )  »  VECT(  1004  )  #.B»  SOi  XR,  MN  ,TSCL,XSCL,SS,SI,  NMX(  2  ) 

NM=  NMX ( N ) 

I F ( ( X P . L  T . X 0 ( 2  »  N ) )  . 0 R . ( X P . G T . X 0 ( N M- 1 t  N ) ) ) G 0  TO  20 


IN  TABLE  N 
END  POINTS 
OF  DATA. 


. 


- 


' 


o  o  o 


C-9 


M=NM-2 
00  10  J  =  2,H 

I  =  J 

I F ( ( X  0 ( J  f  N ) .  L  E  .  X  P  ) .AND. (  X  0  (  J  4- 1 ,  N  )  .GE.XP) )G0  TO  5 
10  CONTINUE 
5  X  M  2 = X  0  (  I  - 1 ,  N ) 

YM2- YO { I - 1 j  N ) 

XM1=X0( I , N ) 

YM1=Y0 (  I  ,  N  J 
XP1=X0( I n ,N) 

YP  1  =  Y0  {  I  H  ,  N  ) 

XP2=X0 { I +2 i N ) 

YP2  =  YC(  H-2, N) 

GO  TO  99 

20  I F { ( XP.LT.XOI 1 , M) ) .OR. (XP.GT .X0( NM, N) ) )G0  TO  30 
I F ( XP  .LT.XOI 2,N) )G0  TO  25 
I=NM-2 
GO  TO  5 

25  1  =  2 

GO  TO  5 


30  WRITE (6, 1002)XP,N 

1002  FORMAT ( '0GJTS1DE  RANGE  IN  TABL4,  XP=',E15.6,'  N=»,I5) 

99  RETURN  / 

END 

SUBROUTINE  SORT  ( I PNTR , YMAX , VMAX , VM I N ) 

FINDS  MAX  DEPTH  AND  WRITES  OUT  T,  X,  V,  0  Y  AT  THIS  PUINT 
FINDS  MAX  DEPTH  AND  VELOCITY  AND  MIN  VELOCITY  FOR  TuTAL  TIME 
THESE  VALUES  CALCULATED  ARE  USED  TO  SET  THE  SCALES  FOR  PLOT 
REAL  MN 

INTEGER  AA,BB,CC 

COMMON  XC(300,2) , Y 0 ( 3 0 0 , 2 )  , T (  25 1 , 2 ) ,X(251,2) , V ( 25 1 , 2 ) ,Y (251,2), 
1M0G( 20,4 ) , VECT( 10J4 ) , B , SO , XR , MN , TSCL , XSCL , SS, SI , NMXI2 ) 

C  COUNT  THE  NO.  OF  VALUES  CALCULATED 

N  =  0 

DO  10  1=1,251 

IF(VECT(4*I-3).LT.0.)G0  TO  15 
N  =  N  4- 1 

10  CONTINUE 

15  AA=4 

DO  20  K=  2 , N 

IF  ( VECT ( 4*K )  .GT.  VECT(AA))  GO  TO  16 
GO  TO  20 

16  AA=4*K 
20  CONTINUE 


BB  =  3 

DO  30  K=  2 , N 

IF  (  V EC  T  {  4  $K.-1  )  .GT.  V  t  C  T ( B  3 ) )  GO  TO  26 
GO  TO  30 
26  BB=4*K-1 
30  CONTINUE 
CC=  3 


DO  40  K=2,N 

IF  (VECT(4*K-1)  .LT.  VECT(CC))  GO  TO  36 
GO  TO  40 
36  CC=4*K-l 
40  CONTINUE 

WR  I  TE  (  6, 2000  )  IP.NTR,  VECTIAA-3),  VECTIAA-2), 
2000  FORMAT (' O'  , 5X , •  STEP  :  ' » I  5 ,10X,  '  T:  ' , F  8 . 3 , 1 0  X , 
1F8.2, 10X, 1  YMAX : ' , F8. 2) 

IF  ( I PNTR  . LE.  1 )  GO  TO  50 

IF  (YMAX  .LT.  VECT  <  A  A ) )  YMAX  =  V ECT ( AA ) 

IF  (VMAX  .LT.  VECT (03))  VMAX=VLCT ( fcB ) 

IF  (VMIN  .GT.  VECf(CC))  VM IN=VECT (CC ) 


VECT ( AA-1 )  ,  VECT(AA) 
'  X: ' ,F3.2»10X, *  V:' 


i 


’ ' 


' 

- 


c-io 


GO  TO  60 
50  YMAX=VECT { AA) 

VMAX=V  ECT ( 03 ) 

VMIN=VCCT(CC) 

60  RETURN 
END 

SUBROUTINE  SCRIB 

C  WRITES  CUT  T  ,  X  ,  V ,  &  Y  FROM  VECT 
REAL  UN 

COMMON  XC(  3  00,  2  )  ,  Y0(  300 ,2),T(25L,2),X(251,2)  , V{251 ,2)  ,Y ( 251 ,2)  , 
1HD0I 20 ,4 ) , VECT ( 1004 ) , b , SO, XR, MN , TSCL , XSCL , SS, SI , NMX! 2) 

C  COUNT  NC.  OF  VALUES  TO  BE  PRINTED  UUT 

N  =  0 

DO  10  1=1,251 

IF( VECT(4*I-3)  .LT.O.  ) GO  TO  15 
N  =  N  +  1 

10  CONTINUE 
15  NX= ( N+ 1 1 ) / 12 
12  =  0 

DO  20  1=1, NX 
11=12+1 
12= 11+11 
I F (  I2.GT.N)  I2  =  N 

WRITE (6,  1000) ( VECT (4*K-3 )  ,K=I1,  12) 

1000  FORMAT (• OT 12F10. 3 ) 

WRITE!€>,  1001)  (VECT!  4*K-2)  ,K=U,  12) 

1001  FORMAT!'  X:',12F10.2) 

WRITE(6,1UG2)(VECT(4*K-1 ) ,K=I1, 12) 

1002  FORMAT!*  V:*,12F10.2) 

WRI TE ( 6, 1003 ) ( VECT (4*K ) ,  K=  I  1 , 12) 

1003  FORMAT!*  Y:*,12F10.2) 

20  CONTINUE 

RETURN 

END 

SUBROUTINE  SHIFT 

C  SHIFTS  THE  INDEX  VALUES  IN  CALCULATED  POINTS 
REAL  MN 

COMMON  XC( 300, 2) , Y0( 300,2) , T ( 25 1 , 2 ) , X ( 25 1 , 2 ) , V ( 25 1 , 2 ) , Y ( 25 1 , 2 ) , 
1  HOG! 20 ,4 ), VECT (  1004 ), B, SO, XR, MN, TSCL , XSCL , SS,  SI , NMX!  2 ) 

DO  10  1=1,251 
T ( I , 1 ) =T ( I ,2) 

X! I ,1)=X(I ,2) 

V! I , 1)=V( I ,2) 

Y( I , 1 )=Y( I , 2) 

IF! T( I , 1 ) .LT.O. )G0  TO  99 
10  CONTINUE 
99  RETURN 
END 

SUBROUTINE  VECTR(N) 

C  STORES  CALCULATED  POINTS 
REAL  MN 

COMMON  XC( 300, 2) , YO! 300 , 2 ) , T ( 25 1 , 2 ) , X! 251 ,2 ) , V ! 251 ,2 ) , Y ( 251 ,2 ) , 
1HDG! 20 ,4 ) , VECT ! 1004 ) , B, SO, XR, Mb , TSCL , XSCL , SS , S I , NMX ( 2 ) 

I PNTR=0 
M=N  + 1 

DO  10  1=1, M 

I PNTR= I PNTR+i 

VECT!  IPNTR)=T ( I,  lJ^TSCL 

IF  (T(  I, L)  .LT.  -0.1)  GO  TO  10 

IPNTR=  IPNTR  +  1 

VECT! IPNTR) =X(I  ,  1)*XSCL 

I PNT  R= I PNTR  +  1 

VECT ( I PNTR) =V( 1,1) 

I PNT  R= I PNT  R+ 1 
VECT! IPNTR)=Y(I ,1) 

10  CONTINUE 


' 


RETURN 

ENO 

SUBROUTINE  U«\ORM(0»N»B,  SO,  SS,  YNORM) 

C  CALCULATES  NORMAL  DEPTH  IN  TRAPEZOIDAL  CHANNEL 
REAL  N 

C  INITIAL  ESTIMATE  OF  YNORM  IS  1.0  FT. 

I  CT  =  1 
YNORM= 1.0 

WO  =  Q  *Q*N*N/ { 2 . 2  2*  SO ) 

S  B= S  Q  R  T ( I . +  S  S  *  S  S  ) *2 .0 
I  Y= YNORM 

A=Y*(B+Y*$S) 

P=B+Y*Se+B+2.*Y*SS 

R=A/P 

W=A*A*R**l .3333 
OW=W-WO 

DY=-3.*A*DW/(W*{ 10.*( B+2.*Y*SS)-4.*R*( SB+2.^SS) ) ) 
YNORM  =YNORM+OY 
TST=ABS( DY/YNORM) 

IF(TST.LE.I.OE-A)GO  TO  99 
ICT=ICT+1 

IF  ( ICT  .LT.  20)  GO  TO  1 
WRITE (6, 1000)  ICT, YNORM 
1000  FORMAT ( 2X , • I CT  =  • , I  5 , 1 OX , * YNURM= • , El 4. 6 ) 

GO  TO  1 
99  RETURN 
ENO 

SUBROUTINE  DC R I T ( 0 , B , S S , YCR I T ) 

C  CALCULATES  CRITICAL  DEPTH  IN  TRAPEZOIDAL  CHANNELS 
C  INITIAL  ESTIMATE  OF  YCR I T  IS  1.0  FT. 

YCR I T= 1 • 0 
W0=Q*Q/32. 17 
1  Y= YCR I T 

A=Y*(9+Y*SS) 

S=B+2.*Y*SS 

W=A*A*A/S 

DW=W-WO 

0  Y  =  -  0  W  *  S  *  S  /  (  A*A*(3  .*S*S-2  .*A*SS  )  ) 

YCR IT=YCRI T+DY 
TST=ABS( DY/YCRIT) 

IF ( TST.LE. 1.0E-6)G0  TO  99 
GO  TO  1 
99  RETURN 
ENO 


' 

■ 

■ 


